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Radio pulsars provide us with some of the most stable clocks in the universe. Nevertheless several 
pulsars exhibit sudden spin-up events, known as glitches. More than forty years after their hrst 
discovery, the exact origin of these phenomena is still open to debate. It is generally thought that 
they an observational manifestation of a superfluid component in the stellar interior and provide an 
insight into the dynamics of matter at extreme densities. In recent years there have been several 
advances on both the theoretical and observational side, that have provided significant steps forward 
in our understanding of neutron star interior dynamics and possible glitch mechanisms. In this article 
we review the main glitch models that have been proposed and discuss our understanding, in the 
light of current observations. 


Radio pulsars are thought to be rotating magnetised neutron stars (NS). The huge moment of inertia (of the 
order of 10^® g cm^) leads to exceptionally stable rotation rates and provides us with some of the most precise 
clocks in the universe. The best timed millisecond pulsars are stable to a precision which rivals that of atomic 
clocks [1]. Nevertheless radio pulsars exhibit several timing irregularities, the most striking of which are the so-called 
glitches. While most objects are observed to steadily spin-down due to the emission of electromagnetic and, possibly, 
gravitational waves, many pulsars show sudden increases in their spin, in some cases followed by an increase in their 
spin-down rate, i.e. glitches. Most pulsars also show slower, long-term, stochastic deviations from a regular spin-down 
law, that are generally classed as ‘timing noise’ and are not the main focus of this review article. 

Soon after the discovery of the first glitches in the Vela nig and Crab HE] pulsars in 1969, several mechanisms 
were suggested to explain these phenomena. Although some initial suggestions featured external mechanisms, such 
as plasma explosions in the magnetosphere [6] or planets around the pulsar [7], the lack of radiative and pulse profile 
changes associated with these events was taken as evidence for an internal origin. The situation is different for rotating 
radio transients (RRATs) and magnetars, which are now known to glitch and do exhibit, amongst other peculiar traits, 
radiative changes associated with these events [EHni- Magnetospheric activity is likely to play a role in these objects, 
as we discuss below. 

There are two main internal glitch mechanisms that have been examined in the literature. The first set of models 
relies on the fact that the outer layers of a neutron star form a crystalline crust that can support stress. As the NS 
spins down, the liquid core adjusts its shape to the rotation rate, while the solid crust maintains the shape appropriate 
for the previous, higher, spin rate. This leads to an increasing amount of strain building up in the crust, which is 
eventually released in the form of a star quake. The quake causes a sudden rearrangement of the moment of inertia 
and ultimately a glitch mM- 

A second set of models is based on the prediction that neutrons in the interiors of mature NSs are superfluid [nms], 
a prediction that has been supported recently by monitoring the cooling of the young NS in the supernova remnant 
Cassiopeia A [nini. At the high densities in the stellar interior (which can exceed nuclear saturation density in 
the core), neutrons find it energetically favourable to pair, thus creating an energy gap below the Fermi surface, that 
must be bridged if the particles are to interact and take part, for example, in the scattering processes that give rise 
to viscosity. Superfluidity thus suppresses the viscosity and lengthens the coupling time-scale between the neutron 
superfluid and the non-superfluid crust. The long recovery time-scales (of order months) observed after the first Vela 
glitch were, in fact, immediately considered to be evidence for a weakly coupled superfluid component in the stellar 
interior [18] . 

Another key property of a superfluid is that its flow is irrotational. The system of neutron pairs is similar to a 
Bose-Einstein condensate, described by a macroscopic wave function, whose momentum is the gradient of a phase. 
Nevertheless we know that the star rotates rapidly and that interactions with the normal component, however weak, 
induce rotation in the neutron superfluid. It is well established from the study of superfluid helium that a superfluid 
rotates by creating an array of quantised vortices that carry the circulation. The angular velocity of an element of 
fluid is proportional to the number of vortices it contains. For it to spin down, some of them must be removed. 
Anderson and Itoh m suggested that interactions between vortices and ions in the NS crust can ‘pin’ the vortices 
and restrict their outward motion. As long as the vortices are pinned, i.e. stay fixed in position, the superfluid does 
not spin down, thereby storing angular momentum, which is periodically released in glitches. 

The vortex model has become the standard picture for pulsar glitches. This is partially due to the crust-quake 
model’s inability to explain the large and frequent glitches now observed in the Vela pulsar |5D], but mainly due to 
the success of the vortex creep model developed by Alpar and collaborators [2T] in explaining post facto the observed 
post-glitch relaxation of the Vela and other pulsars in terms of general parameters of NS structure. Furthermore 
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recent simulations based on this scenario have successfully described several aspects of Vela giant glitches [22] and the 
statistics of the glitching populations as a whole |23ll21]- Nevertheless several issues still need to resolved before the 
vortex picture attains the status of a self-consistent, predictive, falsifiable theory. The trigger for vortex unpinning 
is still unknown and may be due to vortex accumulation in strong pinning areas |25l I26j , vortex domino effects [24] , 
hydrodynamical instabilities [27H25] or quakes 1301 [3I|. Recent calculations have also showed that Bragg scattering 
severely limits the mobility of neutrons in the crust, limiting the amount of angular momentum that can be stored in 
the crust between glitches [32] ■ An analysis of the Vela pulsar reveals that in this case it is difficult to accommodate 
the implied angular momentum exchange during a glitch unless the pulsar has a low mass ^ IMq or part of the core 
is involved in the process [331133] ■ The core of the NS is, in fact, thought to be in a type II superconducting state[TK]. 
Neutron vortices could then interact with proton flux tubes leading to observable consequences [35] . 


I. NEUTRON STAR INTERIORS 

Neutron stars are one of the most extreme and exciting nuclear physics laboratories in the universe. With a mass 
slightly above that of the Sun compressed into a radius of « 10 km, their interiors can easily exceed nuclear saturation 
density, allowing us to probe many-body aspects of the strong interaction at MeV energies, which cannot be studied 
in terrestrial laboratories. The equation of state of matter at such high densities is not known, and several possibilities 
exist. Recent observations of NS masses slightly above 2Mq [36ll37] begin to constrain our understanding of physics 
in this regime. The general relativistic equations of hydrostatic equilibrium, unlike their Newtonian equivalent, lead 
to a maximum mass as a function of radius (or equivalently central density) which is equation of state dependent. In 
figure [^ we plot the mass-radius relation for several equations of state that allow for a maximum mass above 2Mq. 
Softer equations of state (with low maximum mass) are excluded, which constrains the presence of exotica (such as 
hyperons) in the core [55]. 

Furthermore at such high densities the Fermi energy is large compared to the thermal energy so that, although 
mature NSs have an internal temperature ^ 10^ K, they behave essentially as cold objects. In most of the star it is 
energetically favourable for neutrons to be superfluid, while protons in the core are superconducting. In figure [^ we 
show the critical transition temperatures for neutron superfluidity and proton superconductivity for a representative 
sample of models. As can be seen the interior of a mature NS is expected to be mostly superfluid. Note, however, that 
the density dependence of the pairing gaps is such that there will always be regions close to the transition temperature 
in which thermal effects cannot be neglected. 

Superfluidity impacts on the star in many ways. It modifies transport coefficients (e.g. thermal conductivity, 
kinematic viscosity) [391 140] . it allows neutrons and protons to flow independently, and the neutron condensate 
forms an array of quantised vortices (an Abrikosov lattice) which carries the vorticity [15]. We discuss the effects of 
superfluidity in detail in the following sections. 


A. Layered structure 

The outer regions of the star comprise a solid crust, in which ions are arranged in a crystalline lattice which, above 
a density of « 5 x 10^^gcm“^, is immersed in a gas of free (superfluid) neutrons. This region thus behaves as an 
elastic solid and has a finite shear modulus /r, which can be calculated if the equation of state and composition of the 
crust are known [45] : 

4/3 

(Ze)2 g/(cm s2), (1) 

where A„ is the fraction of neutrons outside nuclei, rib is the baryon number density, A and Z the atomic and proton 
number, and e the fundamental charge. For typical compositions and densities [46] this gives a shear modulus of the 
order fi « 10^° g cm“^ s“^. The maximum elastic stress, tm, that the crust can sustain before cracking takes the form 

tm ~ (2) 

where is the breaking strain. By analogy with terrestrial crystals, Smoluchowski and Welch[37] suggested that 
10“® ^ am ^ 10“^, but recent molecular dynamics simulations [35] imply am ~ 0.1. The latter simulations also 
suggest that the system is dominated by pressure and gravity, and the crust does not crack along fault lines, as in 
terrestrial earthquakes, but rather fails as a whole once the maximum strain is exceeded. 

Note that the ions are generally taken to be arranged in a Body Centered Cubic (BCC) lattice, however there is still 
significant uncertainty on the exact composition and structure of the crust [49] . and not only electrons, but also free 
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FIG. 1: LEFT: The mass-radius relation for four representative equations of state: SLy [41], APR [42], GMl [43| and L | 44 |. 
As can be seen the chosen equations of state satisfy the observational constraint of allowing for a maximum mass greater than 
2Mq. RIGHT: The critical transition temperature for a representative set of gap models, taken from [40] (models d,g and 1). 
In the crust one has neutron pairing in the ^S'o channel, with ^P 2 pairing possible at higher densities in the core. We plot also 
the gap for protons in the core. The critical temperature for ®P 2 pairing in the core is constrained to be in the region 
(5 — 8) X 10® K by observations of the cooling of the NS in Cassiopeia A [161117] . The shaded region represents the core. 


neutrons, may partially screen the Coulomb interaction between the nuclear clusters, leading to different, and more 
inhomogeneous, configurations than a BCG lattice [50|[5T]. This can have important consequences for the interaction 
between vortices and ions, as we shall discuss in the following. 

At densities above « 1.6 x lO^^gcm”^ there is a transition to a core fluid of neutrons, protons, electrons and 
possibly muons. The nature of the transition is uncertain. Some equations of state suggest that it is abrupt, while 
others predict a series of phase transitions in which nuclei are no longer spherical but arranged in rods and plates, the 
so called ‘pasta’ phases in which matter behaves more like a liquid crystal than a standard crystal [^. The presence 
of such layers can have a strong observational impact on several phenomena (such as glitches, oscillation modes of 
the star and gravitational wave emissionl [53] not only by lowering the shear modulus but also by smearing out the 
transition region. A perfectly solid crust presents, in fact, a sharp boundary for the interior fluid which reacts to 
changes in rotation by inducing a (dissipative) Ekman flow in a thin layer close to the boundary. The problem is 
complicated by the presence of magnetic fields and by neutron superfluidity but, generally, a smeared out boundary 
leads to weaker dissipation and longer crust-fluid coupling times. 

In the outer core of the star protons are expected to form a type II superconductor [T5] in which the magnetic 
flux is concentrated in flux-tubes. As we discuss below, the strong interaction between superfluid neutron vortices 
and flux tubes can lead to ‘pinning’ also in the core. At even higher densities, above nuclear saturation density, the 
ground state of matter is essentially unknown. At asymptotically high densities the favoured state seems to be a 
condensate of deconfined quarks in the colour-flavour-locked (GFL) phase[M]. If such a phase exists (e.g. in some 
of the heavier stars) it may have interesting consequences; not only is it a superfluid but it also possesses a sizeable 
shear modulus[5S], raising the possibility that ‘core quakes’ contribute to pulsar glitches. 


B. Multifluid hydrodynamics 

In the previous section we purposely glossed over the details of superfluidity in the stellar interior and focussed 
on the constituents at different densities. Let us now consider in detail the consequences of superfluidity. We base 
our analysis on the variational formulation of Andersson & Comer |56j and consider the equations of motion for two 
dynamical degrees of freedom, two macroscopic ‘fluids’: the inviscid neutron condensate and a viscous, charge-neutral 
fluid of protons and electrons. As already mentioned the neutron condensate rotates by forming an Abrikosov lattice 
of quantised vortices, which carry the circulation. Hence a hydrodynamical description only makes sense on length 
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scales larger than, say, the mean intervortex separation, « 4 x 10“^(P/lms)^'^^ cm (see Haskell et al. 2012 for 
an in depth discussion of the length-scales involved). Above this length scale we can average over the velocity field of 
the vortices contained in a fluid element. If the stellar interior rotates uniformly, one can define an angular velocity 
for the condensate, which satisfies the relations: 


KUy 


2[r2n + £n(f^p ~ f^n)] + ~ f^n) 


(3) 


where k, = /i/2m„ = 1.99 x 10“^cm^s“^is the quantum of circulation, with the neutron mass (in what follows 
we assume to„ = nip, the proton mass), Uy is the areal density of vortices, Hn and flp are the angular velocities of 
neutrons and protons respectively, and vj = rsin6* is the cylindrical radius. 

If the individual species are conserved in beta equilibrium, we can write conservation laws for each species. 


9t/0x + Vi(pxO = 0, (4) 

with X = n, p. The Euler equations, on the other hand, take the form: 

{dt + viVj){vf + e^wf") -f Vj(/ix + $) + -b V-’Aj, (5) 

where £x is the relevant entrainment parameter (we recall that EpPp = £nPn), with — uf (y x), 

represents the dissipative part of the stress-energy tensor (see Haskell et al. (2012) [57] for an in-depth discussion of 
dissipation in multifluid systems), px = Mx/wx is the chemical potential and $ is the gravitational potential, which 
obeys the Poisson equation: 


= 47rG^Px- (6) 

X 

The force /f which appears on the right hand side of equation ([^ is the so-called mutual friction force, which exchanges 
momentum between the components and is peculiar to a superfluid system. For straight vortices it takes the form: 


fi = KTlyPnB eijk^iw^y -b KUyPnBcijk^ie'^'-'^niw. 


klr. 




(7) 


where B and B are defined as 


B = 


n 

1 -bP2’ 


and 


b' = 


1 -bP2- 


( 8 ) 


n* is the local angular velocity satisfying Vi = eijk^^x^, and TZ is the dimensionless drag coefficient. 

Mathematically /f has the form of a body force acting on a fluid element, but its correct physical interpretation 
is subtly different. Consider again the most simple case relevant for neutron star interiors, that of a two fluid system 
of superfluid neutrons (denoted n), coupled to a charge neutral fluid (denoted p) of protons and electrons, which are 
locked by the Coulomb interaction on timescales much shorter than the dynamical timescales we are interested in 
[55] . In this situation /f arises from two distinct physical effects. First of all, in the presence of quantised vortices a 
new term appears in the hydrodynamical equations of motion for the neutron condensate. This term takes the same 
form as a (fictitious) Magnus force acting between the vortices and the condensate. 


= unyCijk^^iv^ - vl). (9) 

where vl‘ is the local velocity of a segment of vortex line (averaged over many vortices). Note that the Magnus term 
is only present if vortices are not ‘flowing’ with the condensate. In the absence of additional forces that can lead to 
a relative flow between vortices and condensate (e.g. pinning or dissipation) the Magnus term thus vanishes. Second 
there is a drag force which acts between vortices and the proton-electron fluid. It can be written phenomenologically 

as [55] 


/i = KUyTZivf - A), 


( 10 ) 


where the exact microphysical process that give rise to dissipation will determine the value of the coefficient TZ. In 
the case of straight vortices one can solve for vf by assuming force balance for massless vortices {ff = f^), yielding 
the expression in Q. 

It is important to note that the forces in 0 and ( |I0[ ) are ’average’ forces acting on a small volume around a vortex. 
Given the scales in the problem it is possible to choose, for example, a cylinder such that it’s radius is significantly 
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larger than the vortex core but still small enough to allow us to consider the vortex as a ’line’ in our hydrodynamical 
description. By integrating the momentum flux though the surface of such a cylinder one can obtain the force balance 
condition and derive the expression in 0 by enforcing linear momentum conservation. Strictly speaking, however, a 
vortex cannot be treated as a material object, as its centre is the location of a singularity of the flow, not a massive 
particle or body. The force is not localised on a vortex core, nor does it act equally on all fluid elements within the 
closed surface, as in a rigid body for example. The acceleration of a vortex line is not well defined, so that only the 
average momentum transfer per unit time, and the average force balance condition are well defined. An in 

depth discussion of how the Magnus force arises from the quantisation of vorticity is given, e.g. by Glampedakis et 
al. , Carter et al. [HT] and Sonin [H2] . 

In the core of a neutron star dissipation is mainly from electrons scattering off vortex cores magnetised by 
entrainment [53] (see below) which leads to a drag parameter of the form[59j: 


« 4 X 10"'^ 


/ 6m* 


ml 


1/2 


V0.05/ 


7/6 


lO^^g/cm" 


1/6 


( 11 ) 


where to* = TOp(l — £p) is the effective proton mass, one has iJto* = TOp — to*, and ccp = Pp/p is the proton fraction. 
In the crust, where this mechanism does not act, as protons are not superconducting, the dissipation is mainly via 
phonon interactions |64j or, for larger vortex velocities, via Kelvin waves excited along the vortices themselves |651l66j . 
The values of the dissipation coefficients are, however, quite uncertain. They are likely to be low (10“® ^ ^ 10“^) 

for phonon scattering, while kelvon processes are strongly dissipative (10“^ ^ 7^ ^ 10“^) [53H55] . 

The expression in Q is appropriate for free, straight vortices and can easily be modified to account for the slight 
bending due to the finite rigidity of a vortex (see m for a description). However it is well known from the experimental 
study of He H [55] that a counterflow along the vortex axis destabilises the vortex array to form a turbulent tangle 
(the Glaberson-Donelly instability) [5^. In this case the mutual friction can be described by the phenomenological 
Gorter-Mellink form|70]: 




STT^Pn 

3k 



B 




pn 


( 12 ) 


where « 0.3 and ^2 ~ 1 are phenomenological parameters m- Given that the Reynolds number in a neutron 
star[72] is expected to satisfy Re > 10^, it is likely that instabilities occur, leading to superfluid turbulence and a 
vortex tangle |28j . Even at lower values of Re it is well known from the study of laboratory superfluids that a turbulent 
vortex tangle is likely to develop m- 


C. Vortex pinning 

An important assumption in writing down Q is that the vortices are free to move. However, this is not always so. 
Vortices can ‘pin’ to nuclei in the crust or to magnetic flux tubes in the core of the star, thus preventing the neutron 
fluid from spinning down, as it cannot expel vorticity to keep up with the decelerating protons (to which we assume 
the magnetic field and observed radio pulses are anchored). The pinning force is thus an important ingredient of any 
superfluid glitch model, as it determines the amount of angular momentum that can be stored in the condensate, 
before the Magnus force eventually frees the vortices and gives rise to a glitch. 

Although the pinning force per pinning site is readily estimated, the pinning force per unit length acting on a vortex 
(which is the quantity that balances the Magnus force) is a more complex problem. The first quantity depends only 
on the difference in energy between the state where the nucleus is inside the vortex, and that where it is outside. 
The force per unit length on the other hand depends on several poorly constrained quantities such as the rigidity, 
orientation, and geometry (e.g. tangled versus straight) of the vortex among other things. Several early estimates 
[71H77] considered simple geometries and layers of different pinning strengths. It was also pointed out early on 
by Jones |78| that for an infinite vortex the energy differences between adjacent configurations are negligible once 
averaged over orientations, leading to a vanishing pinning force. Detailed calculations have been carried out by Link 
m, who modelled the crust as a potential due to randomly placed ions, and by Seveso et al. [55] who averaged 
over all orientations of the vortex with respect to a BCG crustal lattice. Both calculations confirm that, for a finite 
vortex rigidity, the force is indeed reduced by the averaging procedure, but is still dynamically significant [SO] El] . 
In figure we show the results of Seveso et al. |80| for the pinning force, also expressed in terms of the maximum 
lag between neutrons and protons that can be built up before the Magnus force breaks the vortices free. As we can 
see the strength of the pinning force spans two orders of magnitude across the inner crust, and has a maximum at 
densities around p « (3 — 7) x 10^^ gm/cm^, depending on the pairing model. 
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FIG. 2: LEFT: The pinning force at different densities obtained by Seveso et al m for their more realistic /3 = 3 model, 
which accounts for correlations in the superfluid pairing gaps. The results are shown for different values of the length scale L 
over which a vortex can be considered rigid. RIGHT: Estimate, for a IAMq NS modelled with the GMl equation of state, of 
the critical difference in rotation rates between superfluid neutrons and the crust AQ = fin — Hp, at which vortices unpin at 
different cylindrical radii (from |80|1. 


We can relate the critical lag in figure (§ to the size of a glitch by a simple conservation of angular momentum 
argument. Let us consider two rigidly rotating components: the crust and all components coupled to it on a short 
time-scale, with moment of inertia /p, rotating with angular velocity flp; and a pinned condensate, with moment of 
inertia rotating at fin, with flp — = A, where A is the critical lag. If all components recover fully after the 

glitch and rotate with common angular velocity fla, one has simply 


/n(flp -b A) -|- Ipilp — {In + /p)fla 


(13) 


The size of a glitch, Ha — Ilp = —/nA/(/n + Ip), thus depends on the lag that builds up and on the fraction of the total 
moment of inertia that corresponds to the pinned superfluid. Note however that post-glitch corotation (i.e. erasure 
of A) implies a correlation between glitch sizes and waiting times, which is not observed in most pulsars (see section 
2). In this sense the above equation should be interpreted as an inequality. Essentially we are saying that in the case 
of strong pinning there is no vortex motion, so no mutual friction force acting between the components. 

Let us examine the problem in more detail by writing the equations in ([^ as [22) : 


tin 

tip 

KUy 


KTlyB 


(Up - II„) _ £n 

(1 - - £p) 1 - £: 

Pn (^p ^n) 


II 


sd 


-B 

Pp (1 ^p) 


-n 

d 


sd 


2[r2n + £n(^p ~ ^n)] + + ^n(^p ~ ^n)] 


(14) 

(15) 

(16) 


where flp and H 


may depend on position, if differential rotation builds up, and Usd is the observed spin-down rate, 
dictated by the external torque T^xt, which is generally assumed to be electromagnetic {T^xt ex flp). Let us first of all 
examine a simplified version of (14)- (dH) in which we neglect entrainment and the external spin-down torque (under 
the assumption that it acts on a much longer timescale than the mutual friction coupling) and let the fluids rotate 
rigidly. The solution is simply: 


flp — rin = Aqc 


1 


I 2Il„B 


(17) 


where Aq = (Up — Hn) is the initial lag in the pinned region, and the coupling time-scale r depends on the ratio 
between the moment of inertia of the proton-electron fluid (and any component locked to it on a faster time-scale) 
and the total moment of inertia I. This is essentially the two-component model proposed by Baym and collaborators 
soon after the observation of the first glitch[TS], in order to explain the slow relaxation observed after the glitch in 
terms of long superfluid coupling time-scales. 
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In neglecting the entrainment, however, we make a serious simplification. In the crust, can be quite large [32j (of 
order £„ « 10). As can be seen from the right-hand side of equation (14), this leads to a spin-down torque acting on 
the pinned neutron condensate, even if we assume ;B « 0 for strong pinning (i.e. zero drag force and mutual friction). 
This strongly reduces the lag that can be built up between glitches and thus the amount of angular momentum that 
can be stored. To see this we can consider the accumulated activity of the Vela pulsar: 




P’ 


(18) 


where the sum is performed over all the glitches and fob is the observation timespan. With the aid of the simple angular 
momentum conservation argument presented above we can relate this quantity to the ratio between the moment of 
inertia of the pinned region and the total that is coupled during a glitch: 

— = —without entrainment (19) 

Ip 

which leads to In/Ip = 0.016, a fraction which can be easily accommodated by the crust. However, recent work 
[551151] has pointed out that in reality one has: 

J 77T-* 

« — A-A — ^ with entrainment ( 20 ) 

Ic lip 

which leads to a larger ratio Is/Ic « 4 — 6 %, using the results of Charnel [55]. This ratio cannot be accommodated by 
the crust, unless the equation of state is soft and the star has mass < IMq. Alternatively it implies that the core is 
involved in the glitch, either by contributing to the pinning or by partially decoupling during and for some time after 
the event. 

As a vortex approaches the threshold for unpinning, the probability rises that it hops to a nearby pinning site 
by thermal activation or quantum tunnelling [82], in what is known as vortex creep. Essentially one can split the 
interaction between a pinning centre and a vortex into a dissipative part, which we model as linear drag, as in the 
previous case, and a non dissipative part, Fpin- The equation for force balance takes the form [55] : 


- Ul) + T^{vf - vf) -b Fp,n = 0, 


( 21 ) 


where Fpin is now the unspecified pinning force. The extra force allows extra freedom to redefine B and B 
independently [51], although note that close to the pinning threshold the linear drag approximation breaks down [75]. 
We discuss creep in detail in section [Vllj However let us note that, in the linear regime which is mostly relevant for 
the short term relaxation of glitches, it can be incorporated in the formalism in (14)-(15), by assuming that only a 
fraction 7 of vortices is free (i.e. are creeping) at a given time, with the others pinned. This leads to a drag force of 
the form = KjnvTZ{vf — vf), which is equivalent to rescaling the coefficient TZ — 'yTZ. We see in the following 
sections that 


Ep (A-Ac) 


■ kT " 


Ac 


( 22 ) 


where A = Hp — is the lag, Ac is the critical lag for unpinning, T is the temperature, k Boltzman’s constant and 
Ep is the energy barrier that has to be overcome for unpinning (that in a first approximation is the pinning energy 

[HI [51]). 

Note that in the above we have assumed that the crust is completely rigid, either due to the rigidity of the lattice 
itself or due to the magnetic field coupling the charged components on short Alven timescales. This is a fairly good 
approximation on longer post-glitch relaxation timescales, but on shorter timescales additional dynamics is likely to 
be present, due to both elastic oscillations of the crust and Ekman pumping of the fluid, even in the presence of 
magnetic fields [55] [55] . 


D. Superconductivity 

Protons are expected to be superconducting in most of the core, with type H superconductivity preferred in the 
outer core, and a transition to type I superconductivity possible in the inner core[60]. For type I superconductivity 
the strength of the mutual friction is still in the range TZ « 10“^ [57]. For type H superconductivity, the magnetic 
field is organised in flux tubes which interact with neutron vortices. The magnetohydrodynamic (MHD) equations 
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for this situation have been derived by Glampedakis et al.|60] and their main feature is that the usual Lorentz force 
is replace by a flux tube tension term. There is then in general an extra force between vortices and flux tubes that 
can lead to vortices pushing flux tubes out [SSJ [53] and, more importantly, creates an energy barrier that can lead to 
pinning in the core. The threshold for unpinning in this case can approach |90j 

/ n \ 1/2 

Aflc ~ 1.5 X 10"^ ( loi^ j lads"^ . (23) 

Once cutting does occur, it is dissipative, with [21] 

/ D \ 1/2 

^ « 2.5 X 10-3 ■ (24) 

It could play a role in several areas, from glitch triggers [55| to their relaxation [3^ [35] and even in limiting gravitational 
wave driven instabilities |91]. 


II. GLITCH PHENOMENOLOGY 


Having discussed the interior physics of NS, let us turn our attention to the phenomenology of pulsar glitches. The 
spin frequency v oi a pulsar is measured by recording the time of arrival (TOA) at the telescope of each pulse and 
fitting a spin-down model obtained by Taylor expanding around a reference time to • 


iy{t) 


Vo + f'o(l - to) + 2^o{t - to) 


(25) 


After subtracting known systematic effects, such as the Earth’s orbital motion and the proper motion of the pulsar. 
Equation (25) is generally a good fit to the data and higher order terms are not usually included. Slow deviations from 
Equation (25) are referred to as ‘timing noise’, while sudden, impulsive, changes in the spin frequency and frequency 
derivative are glitches. 


A. Event statistics 

Glitches were first observed in the Vela and Grab pulsars in 1969 ma (soon after the discovery of pulsars 
themselves |34|) and, while these remain some of the best studied systems, several hundred glitches have now been 
observed in over 100 pulsars [95j . including magnetars and millisecond pulsars. |2 15) The growing volume of events 
allows one to disaggregate the data and highlight the following general trends |33| : 

1. Glitch sizes and waiting times. Glitch sizes span several orders of magnitude, with the smallest glitches of 
the order of Avjv « IQ-^^ and the largest of the order of Avjv « IQ-^ (see e.g. Espinoza et al. [95]). Melatos 
et al. |96j showed that the size distributions are consistent with power laws, with the index varying from pulsar 
to pulsar, and the waiting time distributions are consistent with exponentials, as one would expect from a scale- 
invariant or self-organized critical process, such as earthquakes or vortex avalanches m- There are, however, 
significant exceptions to this trend: in Vela, J0537-6910 and J1341-6220, glitches recur quasiperiodically [361138] 
and exhibit a narrower spread in size. This suggests the presence of an additional mechanism that sets some 
natural scale for the process and may set the maximum scale for glitches, although a rapidly forced critical 
system may also exhibit quasi-periodic avalanches, as we discuss in the following section. Very interestingly, 
a recent analysis of the lower end of the Crab’s size spectrum, suggests the existence of a minimum size of 
Av jv « 10“®, which is significantly above the smallest resolvable event [33]. The lower end of the distribution, 
however, is contaminated by a different population of timing noise events (possibly of magnetospheric origin) 
and further studies are required to identify a lower limit securely. At the upper end of the spectrum, across the 
entire population, Espinoza et al. [95] have identified a sub-class of ‘giant glitchers’ that exhibit large glitches 
of similar sizes, Av « 10/i Hz. 

2. Reservoir effect. Do larger glitches occur after longer waiting times, i.e. does the star contain an angular 
momentum ‘reservoir’ which is replenished before and is fully emptied during a glitch? Several models (e.g., 
star quakes) predict a size-waiting-time correlation, whereas others (e.g., self-organized critical systems) do not. 
Strikingly, such a correlation exists in the pulsar J0537-6910. However, the correlation does not exist in the rest 
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FIG. 3: LEFT: The correlation between glitch sizes and waiting times to the next glitch. RIGHT: The correlation between 
glitch sizes and waiting times from the previous glitch. The lines indicate the best linear fit and the correlation coefficients are 
shown in table [^.The only case in which there is a correlation is for the pulsar J0537-6910, in which there is a correlation 
between the size of a glitch and the waiting time to the next, suggesting that the size of the glitch is random, but sets the time 
it takes to build up enough stress to reach the threshold again, as in the crust quake model. 


TABLE 1: The correlation coefficients for the log-log hts in hgure for correlations between glitch size and waiting time to 
the next glitch (case 1), and between size and waiting time front eh previous glitch (case 2). As can be seen on the pulsar 
J0537-6910 shows a significant correlation in case 1. 


Pulsar 

P(l) 

P(2) 

J0537-6910 

-0.018 0.94 

J0835-4510 (Vela) 

0.72 

0.017 

J0534-b2200 (Grab) 0.42 

0.032 

J1341-6220 

-0.062 0.63 

J1740-3015 

-0.083 0.24 


of the population; in general, there appears to be no ‘reservoir’ effect [35]• In figure]^ we plot glitch size verses 
waiting time to both the previous and next glitch. Only J0537-6910 shows a correlation between the size of a 
glitch and the waiting time to the next, which can be interpreted as the time to reach a critical stress being 
proportional to the size of the previous glitch |100j . Melatos et al. [96] argued that the correlation observed 
in J0537-69I0 is, in fact, a by-product of the observed quasi periodicity and the fact that large glitches, of 
approximately the same size, dominate the size distribution. 

3. Age and glitching activity. Approximately 10% of the pulsar population exhibits glitches, with seven pulsars 
having glitched more than 10 times . In general, middle-aged pulsars glitch the most; activity decays with age, 
and older pulsars have smaller glitches m- There also appears to be a correlation between glitching activity 
and the spin down rate [35]. The Root Mean Square (RMS) variations in the residuals associated with timing 
noise also correlate positively with the spin-down rate [Ml IM]- Electromagnetic spin down drives glitch 
activity in most models, so it is natural that older pulsars (which spin down slower) are less likely to glitch (let 
alone have large glitches) during the observations [35] . 

For pulsars that have glitched often another interesting quantity is the amount of the spin down that has been 
reversed by glitches, the activity parameter: 




(26) 


As already discussed in the previous section this can be linked to the moment of inertia of the region initiating 
the glitch, from relations (191 and (20). In table we show the activity parameter, and constraints on the 
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TABLE II: We list, from reference [34], the characteristic spin down age Tc and the activity parameter A for several pulsars, 
as well as the f ract ion of the total moment of inertia that is transferred during the glitch (the superfiuid ’reservoir’), obtained 
from equation (191. Note that if one accounts for strong entrainment and uses equation (201 to evaluate the ration of moments 
of inertia, the results are generally a factor 4-6 larger. 


Pulsar 

Tc (kyr) 

A (10- 

-Vd) IJI (%) 

(no entrainment) 

J0537-6910 

4.93 

2.40 

0.9 

B0833-45 (Vela) 11.3 

1.91 

1.6 

J0631-tl036 

43.6 

0.47 

1.5 

B1338-62 

12.1 

1.31 

1.2 

B1737-30 

20.6 

0.79 

1.2 

B1757-24 

15.5 

1.35 

1.5 

B1758-23 

58.4 

0.24 

1.0 

B1800-21 

15.8 

1.57 

1.8 

B1823-13 

21.5 

0.78 

1.2 

B1930+22 

38.8 

0.95 

2.7 

J2229-b6114 

10.5 

0.63 

0.5 


moment of inertia of the glitching region, taken from Andersson et al. [Mj. As we can see in general one needs 
at least a few percent of the total moment of inertia of the star to be transferred during glitches, possibly up to 
nearly 10% if we account for the large neutron effective masses calculated by Charnel [32| . 

4. Radiative and pulse profile changes. The general lack of correlation between glitches and radiative changes 
was taken early on as evidence for the internal origin of these phenomena. Nevertheless there are now some 
high field pulsars for which radiative changes, such as changes in the pulse profile, are observed in coincidence 
with glitches IHHII] as is the case for magnetars |104j . It is thus clear that the magnetosphere plays a role in 
these systems, although how this could be, when the dynamical time-scale of the magnetosphere (~ spin period) 
is much shorter than the glitch waiting times, remains a mystery. In the case of magnetars, it is also thought 
that sudden changes in the magnetosphere are behind the observed phenomena of ‘anti-glitches’, i.e., sudden, 
negative Av. Correlated changes in the pulse profile and spin-down rate are also associated with timing noise 
events in several pulsars [8]. 

A recent analysis of the Crab’s spin-down rate and glitching activity has also claimed a link between an 11 year 
period of increased glitching activity, and changes in the breaking index in the same period, which would be 
linked to a change in the spin-down torque |105j . 


B. Post-glitch recovery 


The rise time of a glitch cannot be resolved at present. The best upper limit is « 30 s for the Vela pulsar if 
one neglects the two ‘slow glitches’ observed in the Crab, which are thought to have been a different kind of event 
[10611107] . After the unresolved rise, there is a change in the spin-down rate, part of which is permanent, and part of 
which decays over days. It is customary to fit these relaxation events with a function of the form: 


N 

v{t) = Vsd{t) + Avp + Avpt + ^ , 

i=l 


(27) 


where I'sdit) is the pre-glitch spin-down law in (25), Avp is the ‘permanent’ step in frequency. Ah' is the ‘permanent’ 
step in the spin-down rate, and Aui are the amplitudes of the decaying parts of the frequency step, with decay time- 
scales Ti. The number of exponentials N that is fitted varies from pulsar to pulsar and glitch to glitch and depends 
on how long after the glitch one has observations and on the sampling rate. For some of the better observed Vela 
glitches one needs to fit four exponentially decaying terms [108] . The total size of the glitch, Auq, is 


N 

Avq = Avp 

2=1 


(28) 
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The relaxation can vary significantly from pulsar to pulsar, and even from glitch to glitch in the same object. The 
latter property is especially surprising: glitches are small-amplitude events with ^vcjv 1, so one naively expects 
the recovery to be a linear process, whose time-scale is independent of its amplitude and depends on constitutive 
coefficients (e.g., viscosity and mutual friction) whose values do not change appreciably during inter-glitch intervals of 
years [SS]. In some objects the recovery is well approximated by simple steps in frequency and frequency derivative, 
but often this is due to infrequent observations that make it impossible to fit for short-term exponentials. Other 
glitches appear to be simple steps in frequency alone. Where is actually measured, it is sometimes unclear 
whether one is measuring the enhanced spin-down rate during a long recovery, i.e., Az> « or a genuine 

permanent torque change Aj/j,; it is a priority to clarify this issue observationally. Often different kinds of relaxation 
are observed in different glitches of the same object [Ml [IDS]. Furthermore ‘microglitches’ of small amplitude can have 
all kinds of signature, with negative and positive steps in v and v mniiiii]. Recently a population of such events 
has been uncovered in the Crab pulsar [99] and appears to be separate from the main glitch population and of likely 
magnetospheric origin (i.e., to be, in some sense, part of the timing noise). 

The best resolved glitches are those of the Vela and Crab pulsar, for which we have the best coverage and largest 
body of data. The recoveries in these two NSs are quite different. For the more recent glitches of the Vela pulsar it is 
been possible to measure large increases in the spin-down rate, Ai>/i> % 1, which cannot be accounted for simply by 
assuming that part of the crust is decoupled but require changes in the internal torque ina. Moreover, if we define 
a healing parameter: 


q = E 


Ai/j 

Az/g’ 


(29) 


one finds Q > 0.8 for the Crab |107j . and Q < 0.2 for Vela [113] . In general Q appears to correlate with v and young 
pulsars like the Crab exhibit near complete recovery and older pulsar like Vela only a very partial recovery [114] . The 
measurements in Vela are complicated by the fact that a significant fraction of the glitch has not recovered before 
the next one occurs. It is interesting to note that in the Crab the persistent offset in the spin-down rate (Az>p)leads 
to the star ’overshooting’ and spinning slower than it would have had the glitch not occurred (an effect that may be 
masked in Vela by the occurrence of the next glitch). Link et al. [115] argued that this is not possible in standard 
superfluid creep or crust quake models, unless one allows for a time-dependent torque that changes during the glitch 
(i.e., that is not purely a function of the lag between the crust and the superfluid) or if there are changes in the 
external torque, e.g., due to magnetic field rearrangement. Van Eysden & Melatos (2010) |H5| provided a different 
perspective, showing that the overshoot emerges naturally in a two-component superfluid model of glitch recovery 
featuring Ekman circulation and mutual friction, provided that the Ekman and mutual friction time-scales lie in a 
specified range. The van Eysden-Melatos model reproduces spin-down in laboratory helium experiments to better 
than 1% [116] . 

The recoveries of magnetar glitches are also remarkable in many ways. Many exhibit large increases in the spin- 
down rate {Avjv « 0.5) which persist over weeks to months [11711118] . Haskell & Antonopoulou [119] suggested that 
this is a consequence of the long rotational period of these systems and that the physics governing the recoveries in 
magnetars is the same as in radio pulsars, although the glitch trigger may not be, with magnetospheric effects likely 
to play a role. It has, in fact, been suggested that torque variations associated with magnetic field reconfiguration in 
the magnetosphere may be at the origin of the peculiar ’anti-glitch’ that has recently been observed in the magnetar 
IE 2259-1-586 by Archibald et al. [12011121] . In this event there appears to be a sudden spin down (rather than up as 
in conventional glitches) of the star, followed by another timing event, which is harder to interpret, but is most likely 
a second anti-glitch [122] . 

Upper limits on the temperature increase of the Vela pulsar have also been obtained from Chandra X-ray data, 
limiting the increase in temperature to < 0.2% 35 days after the glitch [123] and < 0.7% 361 days after the glitch 

[mj. 


C. Timing noise 


Timing noise is a general term to define stochastic deviations from equation (25). It manifests as noise in phase, 


frequency, or frequency derivative in different objects [12511127] . In general the effects can be due to various processes 
includine precession [12811129| . Tkachenko modes m, superfluid vortex creep |M]) superfluid turbulence [103] . or 
changes in the magnetic field [1011 IllOl 1131] . Timing noise is of interest for the current discussion as it has been 
suggested that glitches below the observational threshold may be the cause of timing noise [Ml ina and of the 
measurement of anomalous breaking indices [13311134] . Janessen & Stappers [135] in fact showed that small glitches 
can account for part of the timing noise in PSR B1951-132, but not all of it, while Hobbs et al. m claimed that 
for young pulsars with Tc < 10® yr timing noise is dominated by glitch recoveries. However recent measurements 
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have highlighted correlations between pulse shape and timing noise [5] , while surveys [95119911136] have uncovered a 
population of small, impulsive events that can be fitted as glitches as long as one allows for negative and positive 
Aiy' and At'. This population appears to be distinct from the larger glitch population and is likely of magnetospheric 
origin. If this is the case, glitch statistics are confusion-noise limited, as one reaches the level of sensitivity necessary 
to uncover timing noise 


Triggers 

The rest of the review deals with the theories put forward to explain the observations summarized in Sec. 2. We 
divide the theories roughly into two groups: trigger mechanisms, which explain how the glitch originates in the first 
place and predict glitch sizes, waiting times, and their distributions; and relaxation mechanisms, which explain the 
shape and time-scales of the recovery. This distinction is convenient but somewhat artificial, as some theories address 
both issues, as discussed below. 


III. CRUST QUAKES 


Soon after the first observations of glitches in the Vela and Crab pulsars, crust quakes were suggested as possible 
causes for these phenomena [HIS]. The fact that terrestrial earthquakes follow a power-law energy distribution 
makes star quakes an attractive explanation for the glitch size distribution observed. The crust forms a rigid lattice, 
which crystallized at a previous, faster, spin rate. As the star spins down, the crust maintains this shape, but it is no 
longer the equilibrium shape corresponding to the current angular velocity. Stress thus builds up. To quantify this 
process let us define the oblateness as m- 

e = {I - Inr )/Inr, (30) 


where / is the actual moment of inertia of the star, and Inr is the moment of inertia of the star in the absence of 
rotation. The equilibrium value of e is determined by the competition between gravity, which tends to make the star 
spherical, and rotation and elasticity, which tend to make it oblate. The gravitational energy takes the form: 

Eg = E^'" + Ae\ (31) 

where E^^ is the energy of the non-rotating configuration, and the coefficient A is determined by the stellar structure; 
for an incompressible sphere one has A = (3/25)(GM^/ii). The deviations are quadratic in e as a sphere is the 
minimum energy configuration. The elastic energy can also be expanded around the unstrained state, that corresponds 
to a reference value cq. In this case we have 

Eel = Bie - eof, (32) 

with B = (57/50^) (47ri?^/3), for a self gravitating sphere of constant shear modulus fi, and, more generally, B A. 
The rotational energy is 

Er = L^/2I , with L = m, (33) 


where L is the total angular momentum. Upon minimizing the total energy Eg + E^i + E^ while keeping L constant, 
one finds 


Inr^^ , B 
A{A + B) A + b'^°' 


(34) 


The reference oblateness is 


eo 


AA 


(35) 


with Qq the rotation rate corresponding to the relaxed state (i.e. the rotation rate at which the crust would be 
completely un-stressed). The mean stress in the crust is f = ^(cq — e) When a critical strain ttc = A/m is reached, the 
crust readjusts to a new reference shape (a crust quake) with a sudden negative change Acq. This shifts e (the actual 
ellipticity of the star) by an amount: 


B 

A + B 


Aeo, 


Ae = 


(36) 
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which leads to Ae = ^ = where we use equation (30) to relate the change in oblateness to a change (decrease) 
in moment of inertia and thus to the observed frequency jump All. 

Given the size of a glitch, the quake model predicts the size of the next event. The stress relieved during a quake 
is: 


At = /r(Aeo - Ae) = M 

ID 


while the rate at which strain builds up is 


t = —fii = 




/ 02 


4^1 Tsd 


where we define the spin-down age by Tsd = —fl/n. The waiting time to the next glitch is thus: 


Sta 


t 


4:A^ Tsd All 

B fl ■ 


(37) 


(38) 


(39) 


Note that this assumes semi-complete release, which never happens in terrestrial earthquakes and Self Organised 
Critical Systems (SOCs) in general [Ml I138j . Furthermore the change in moment of inertia produces a permanent 
change in the spin-down rate (assuming no change in the external torque) 


Afi 

n 


A/ 


Afi 

IT' 


(40) 


We can now compare the predictions to observations. The main difference from one star to another could, of course, 
be the factors A and B, which depend on mass and radius. Nevertheless, for typical masses and radii (M « 1.4Mq 
and i? « 12 km) the large glitches in the Vela pulsar would require almost all the stress accumulated in the crust from 
birth (assuming a period at birth of 1 ms) to be released during the glitch and furthermore would require such events 
to be rare, with waiting times of tens of millions of years. This is clearly not in accordance with what is observed; 
Vela has large, frequent glitches (Afi/fi « 10“® approximately every 2-3 years). The situation is less severe for the 
Crab, but the core quake model would still require the Crab to have a low mass (< IMq) [13711139j . Equation (40) 
also predicts smaller permanent steps in the spin-down rate than are observed in Vela, and a correlation between Afi 
and An that is not generally observed in the pulsar population |96j . There is also no observed correlation between 
the size of a glitch and the waiting time until the next, with the notable exception of PSR J0537-6910 (one of the 
quasi-periodic glitchers), where a linear correlation has been claimed [100] . 

The main ingredient missing from the quake picture is neutron superfluidity. Superfluidity plays an important role 
in building up stress in the crust, leading to quakes |140H142] . This can be seen if we consider the average pinning 
force per unit volume acting on the lattice. 


= nyfi{w,z). 


(41) 


where ft, the pinning force per unit length that opposes the Magnus force, is directed along the cylindrical radius zu 
for a straight vortex. As fi depends on density, it varies along a vortex, leading to 

(42) 

Since Ff shears the lattice it cannot be balanced by gravity, pressure or centrifugal forces, which are all expressed 
as gradients, but must be balanced by elasticity, building up a strain. The maximum lag that can be accumulated 
before the crust breaks can be estimated as [140j : 


Afi 


_ ^lGM 

'^''pmi 


(43) 


with CTc the critical breaking strain, p the pressure and fi the angular velocity of the star. If the crust is weak 
(be « 10“^) one gets reasonable values of the lag, which could be reached in ~ 3 yr in Vela. However, the more 
realistic breaking strain, CTc ~ 0-l> obtained with molecular dynamics simulations |48j implies that vortices unpin well 
before the critical strain is reached, ruling out crust quakes as a trigger mechanism for large glitches (but not for 
smaller glitches and in younger pulsars such as J05370—6910 or the Crab) [14311144] . More detailed modelling of the 
breaking of the crust needs to be undertaken to confirm the nature of the fracturing process [ISl I145| . Carter [146] 
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showed that stress can build up even in the absence of vortex pinning, as long as the viscosity between the neutron 
superfluid and the crust is low. In this case, if the neutrons rotate differentially, the crust must be strained to balance 
the Magnus force, and the mechanism proceeds as described above, generating star quakes once the breaking strain 
is exceeded. 

Ruderman mz] has suggested that stress may also build up as superconducting flux tubes are pushed out by 
vortices during the spin down of the pulsar. The footpoints of the flux tubes are attached to the crust and shear it, 
leading to a quake. During the quake, platelets may break off, carrying the frozen-in magnetic field and migrating 
towards the equator, where they are subducted. This reduces the magnetic dipole moment and hence the external 
torque, accounting for a permanent change in v. A detailed analysis, however, suggests that vortex motion would 
mainly wind up the field in the azimuthal direction |148) . The model suffers from the same problems described above, 
namely that the crust would not break before vortices unpin (even if pinning is to flux tubes), and the crust cannot 
support cracks such as those needed for platelet migration [351 1149j . Furthermore the model predicts that pulsars 
with ^ < 1 Hz should not glitch (as essentially all the magnetic flux has already been pushed out, if they were born 
spinning fast), which is not the case. 

Another possibility is that the crust is indeed too weak to be the seat of glitch activity, but that the star harbours a 
solid core which experiences ‘core quakes’ |150j . Recent calculations suggest that the core may, in fact, be composed 
of a quark condensate in the colour-flavour-locked (CFL) phase |151j . with a shear modulus ^ « 10^^ — 10^“^ erg/cm^ 
[152] . The greater shear modulus and moment of inertia can lead to larger glitches and lower waiting times. In fact it 
has been noted by Alpar [144] that one may expect larger glitches than are observed in this case. More work is needed 
to understand how a solid core fails upon reaching the critical stress, if weaker sections exist, the overall nature of the 
phase transition, and the predicted statistics. 


IV. VORTEX PINNING AND THE ’SNOWPLOW’ MODEL 


The standard view of pulsar glitches is that there is a superfluid component (usually assumed to be in the crust) 
whose vortices are pinned. Then, the superfluid cannot lose vorticity and spin down, lagging behind the normal fluid, 
which spins down electromagnetically along with the solid crust [19] . An important ingredient in this scenario is the 
maximum pinning force that the lattice in the crust (or flux tube array in the core) exerts on the vortices. Let us 
expand the discussion on the pinning force from section (IC) . Recent estimates that account for finite rigidity and 
different orientations of the lattice allow a sizeable lag to build up, large enough to explain Vela glitches [80]. Let us 
now analyse this scenario in more detail. 


In his seminal work in 1970 Alpar m estimated the maximum pinning force acting on a nuclear cluster as: 


\Fp 


no 


A‘^{nG) A2(no) 


Epinc) °~Ep{no) 


V 

X 


(44) 


where A is the supefluid gap, Ep the Fermi energy of the neutrons, V is the volume of a nuclear cluster, no is the 
neutron number density inside the clusters, while no is the number density outside (in the free state), and x is the 
scale of the interaction. One has x = min(i?vys, ^), where Rws is the radius of the Wigner-Seitz cell, ^ = Ep/{kpA) 
is the vortex coherence length (the radius of the vortex core), and kp the Fermi momentum of the free neutrons. 
Essentially this is the difference between the energy cost of a vortex core of normal matter being placed at a density 
no rather then no- If this is energetically favourable, with A"^(nc)/Ep{nG) > noA‘^{no)/Epirig), vortices pin to 
clusters. To obtain the pinning force per unit length, to compare to the Magnus force, one simply divides by the 
lattice spacing a = 2R\ys- This estimate, and subsequent analysis such as that of Epstein and Baym [153] who 
included the contribution of the kinetic energy of the superfluid, leads to rather large forces, of the order of fp = 10^® 
dynes cm“^, corresponding to a critical lag of Ail^ « 10 rad s“^. The critical lag is much larger than the lag that 
builds up between glitches (All « 0.01 rad s“^) for Vela, easily accommodating the angular momentum transfer 
needed to explain the glitch. It is also consistent with the lack of a reservoir effect in most pulsars, i.e. with the fact 
that glitches do not appear to release all the angular momentum that has been stored since the previous event (see 
section |n] for a detailed discussion). It requires a trigger mechanism that would release the stored angular momentum 
well before the lag corresponding to the maximum of the pinning force is exceeded. Several suggestions have been 
made, including quakes freeing vortices [15411155j (although, as discussed in the previous section, it is unlikely that 
the crust yields before the vortices unpin) and pinning inhomogeneities leading to regions of vortex over-density and a 
locally perturbed Magnus force [551 1142j . Gross-Pitaevskii simulations do, in fact, show that over-densities of vortices 
are key drivers of vortex unpinning 

The force obtained in equation (44) is large mainly due to the simplified geometry, which assumes that the force 


per nuclear cluster acts over the whole vortex, essentially assuming that the vortex is aligned with one of the principal 
axis of the crustal crystal. Jones [75] pointed out that if one averages over all possible orientations of the crystal 
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(effectively treating the crustal lattice as a random potential, which would also be the appropriate treatment if the 
crust is an amorphous solid), the energy difference between the pinned and unpinned states depends on the fractional 
difference in the number of nuclei threaded by the vortex in the two states, leading to a vanishing pinning force in 
the limit of infinitely long vortices. 

The length scale over which a vortex can bend to accommodate pinning sites has been estimated by Link and 
Cutler |157j as I « (10^ — \Q‘^)Rws- Seveso et al. [80] obtain a slightly larger value of I « (10^ — 10^)i?vvs and a 
different functional dependance on the tension and pinning energy, in agreement with the calculations of Hirasawa 
and Shibazaki [158] . Given the length scale over which a vortex can be taken to be rigid, one can calculate the pinning 
force per unit length. This has been done by Link Hi mi, who simulated vortex motion in a random potential and 
obtained a pinning force per unit length fp of the form 


f ^ / AEp \ 

^ Rwsi K^R-wsTv) 


(45) 


where Ep is the pinning energy and T„ is the self energy (tension) of a vortex. For typical numbers this gives forces 
of the order fp « 10^® dynes cm“^. 

Recently a similar value was obtained by Seveso et al. [80] . who calculated the pinning force per unit length of a 
vortex, using the superfluid gaps of Donati & Pizzochero |159j . In this case the the force is obtained by integrating 
over all vortex orientations the quantity 

= (46) 

where An{9, (f>) is the difference in number of intercepted nuclear clusters, as a function of the angles (0, (j)) which 
represent the orientation of the vortex with respect to a reference position. Averaging over all possible orientation 
reduces the force considerably, leading to forces in the range fp ~ 10^® dynes cm“^, and critical lags of Afl ~ 0.01 rad 
s“^. The exact value depends on the region of the crust considered and weakly on the equation of state and stellar 
model. Figure shows an example of a realistic density dependence of the lag, which clearly has a maximum in 
the strongest pinning regions. These more realistic values are comparable to the lag that Vela can build up between 
glitches, and the appearance of a maximum in the equatorial region, allows for a new picture of glitches for pulsars 
which exhibit a reservoir effect, the so-called ‘snowplow’ model |5S]. 

In the snowplow model vortices unpin in the stellar interior, as the threshold for unpinning is exceeded, and are 
assumed to eventually repin in the stronger pinning regions at the equator, gradually moving outward as a larger lag 
builds up. The model assumes that a glitch occurs once the critical lag corresponding to the maximum of the pinning 
force is reached. At this point the vortices are distributed in a thin sheet close to the region in which the pinning 
force peaks, containing Ny = (2Tr / Aft max vortices, where Vmax is the cylindrical radius at which the maximum 

A^max of the critical lag is located. Once the maximum is exceeded vortices are free to move out an annihilate at 
the edge of the neutron drip point Rnd, leading to an exchange of angular momentum: 

rRnd _ 

ALgi = 2KNy / dxx / dzpa{y x'^ + z^)., (47) 

Jo 


where Pn{i") is the radial density profile of the superfluid neutrons. To obtain the physical observables of the glitch 
one must make assumptions about the superfluid fraction of the moment of inertia = In/I and the fraction of 
total moment of inertia that remains coupled to the ’normal’ component of the crust during a glitch, denoted Ygi. 
The step in frequency is then given bv[26] 


Ai' = 


ALgi 1 

27rl 1-Q^(l-Ygi)’ 


(48) 


and the step in frequency derivative, immediately after the glitch, by 


^ ^ Qn{l-Ygi) 

V 1-Q„(l-Vg0' 


(49) 


The waiting time between glitches is simply given by the time to build up the maximum lag again after a glitch: 




n 


(50) 


By fitting the above expression to the average waiting time of Vela glitches (which we remind the reader is a quasi- 
periodic glitcher) one can then fix the amplitude of the pinning force (up to uncertainties in the vortex tension) and 
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then use (491 to estimate Ygi , having calculated Qn from a set of representative equations of state |160) . By comparing 
the predictions to Vela glitches one can then estimate that less than 10% of the moment of inertia must be decoupled 
during the glitch. By comparing equation (49) to the short term increase in spin-down rate after the 2000 and 2004 


glitches one can constrain the gap model used for the calculations |160j . A comparison of equation (501 to the average 
waiting times for giant glitches in the sample of Espinoza et al. [SS] reveals good agreement also for these systems 


Despite these encouraging results there are some issues with this model. The first is that the equations in ( 48]|^ ) 
do not account for entrainment in the crust |32l I33| , which strongly reduces the lag, and hence the amount of angular 
momentum that can be exchanged. This effect can be parametrised in the model by allowing for larger values of 
Ygi. However these are not easy to reconcile with the constraint in (491 and the activity of Vela [161] . Haskell et al. 
[22j simulated the hydrodynamical response of the star to a glitch triggered by the snowplow mechanism and showed 
that, for some hours after the glitch, the main contribution to the increased spin down rate comes from the mutual 
friction in the outer core, not from the external torque. The constraint in (49) is thus not relevant for the observed 
short-term increase in i> of the Vela pulsar. Haskell et al [93] also show that to fit the post-glitch relaxation in Vela 
one cannot allow for part of the vorticity in the core to remain pinned to flux tubes after the glitch, as this would 
lead to the crust relaxing slower than observed. 

A more important issue is that the snowplow model predicts glitches to be periodic, contrary to what is observed in 
most pulsars. So while the model can explain the glitching behaviour of Vela and the other quasi-periodic glitchers, it 
fails to account for other systems in which the size distribution is well described by a power law and the waiting time 
distribution by an exponential |96j . In general it is thus expected that, while the snowplow model can set the maximum 
size of glitches in a system, other mechanisms must be active to account for the smaller, more frequent glitches. Quakes 
may be responsible, as already discussed, but an intriguing alternative is that vortices move collectively in avalanches, 
as we now describe. 


V. SCALE INVARIANCE AND VORTEX AVALANCHES 

The apparent power law statistics of pulsar glitch sizes, spanning up to four decades in individual pulsars, point to 
a scale invariant process. One possibility that is especially pertinent to a slowly driven vortex array is Self Organised 
Criticality (SOC). Self organised criticality was first invoked by Bak |162[ 1163] to explain how simple physical rules 
can lead to complex behaviour in sub-critical systems (i.e. systems which are marginally stable and close to an 
instability threshold). Examples of SOC can be found in several fields of physics, ranging from earthquake and forest 
fire modelling to superconductors |92| • Warszawski and Melatos |23l [24l 1164] have shown that pulsar glitches may by 
a manifestation of SOC, which would explain many aspects of the phenomenon and describe how collective motion 
of vortices can trigger glitches of varying size. In particular two different mechanisms have been considered as the 
source of scale-invariant dynamics in pulsars: vortex avalanches [531 [53] and coherent noise processes |165| . Let us 
analyse both these mechanisms in more detail. 


A. SOC and Nearest-Neighbour Avalanches 

A system in an SOC state is generally composed of many discrete elements, whose motion is dominated by nearest 
neighbour interactions rather than by global forces. Each element moves when the local force exceeds a threshold, 
thus allowing for stress to accumulate in certain regions and to relax rapidly in others. The system evolves through 
a series of metastable states, driven by an external driver which operates on a much slower time-scale than that 
associated with the local interactions. Over long time-scales, the system self-adjusts to fluctuate around a ‘self-tuned’ 
configuration which is at the threshold of relaxation, e.g., a sandpile tends to the slope at which sand grains start to 
topple [138] . In the latter sense, the system is ‘critical.’ 

This definition resembles quite closely our understanding of superfluid vortex dynamics in a pulsar crust. Elec¬ 
tromagnetic spin down drives the system on long timescales. Vortices unpin when the Magnus force exceeds a local 
threshold, and interact with each other locally on a fast time-scale. The critical state here is exactly analogous to 
the Bean state in a type H superconductor |166| : a vortex array which is just irregular enough such that each vortex 
is on the threshold of unpinning. (Note that this is not a regular Abrikosov lattice, where the Magnus force is zero 
everywhere). Alpar suggested early on [551 116711168] that there may be capacitive and depleted regions of vorticity 
in the crust, in which steep gradients in the vortex density produce an increased Magnus force. Random unpinning 
may then lead to a vortex avalanche via knock-on effects. In fact, this is exactly what is observed in SOC systems. In 
general transitions from one metastable state to the next occur via avalanches that have no preferred scale and lead 
to distributions of sizes and lifetimes (durations of the avalanches) that are power laws. 
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To investigate this scenario quantitatively Warszawski and Melatos |164) constructed a two-dimensional cellular 
automaton, containing N'^ cells covering an equatorial section of a star of radius R. To make the problem tractable, 
vortices are grouped in bundles, which are assumed to be the discrete elements of the simulation. In general there 
are = Ny/{Bi,N^) elements, where Bb is the average number of bundles per cell, and Ny is the total number of 
vortices in the star, obtained from equation ([^. Note that only e7V„ vortices are assumed to be pinned; the rest keep 
up with the externally induced spin down of the crust. This feature of the model is unrealistic, at least in the inner 
crust, where pinning sites are separated by 10^ fm. One then defines a maximum pinning force, which gives a 
threshold Magnus force for unpinning Fth, assumed to be uniform. For each cell one then compares this threshold to 
the Magnus force acting on the bundles, which is calculated as 

= (51) 

where v'" is the vortex velocity (equal to the crust velocity, as the vortices are pinned), and the velocity of the 
superfluid neutrons is calculated from 


^pinned "^^unpinned ^NN > 


(52) 


where the first two terms are due to the pinned and unpinned vortices contained within a circle of radius rcen, the 
position of the cell, and calculated from the Feynman relation in equation (|^. The last contribution is that of the 
eight nearest neighbouring cells. The system is evolved in time according to the following rules: 


1. The grid is initialised by placing vortex bundles at random. 


2. The Magnus force is calculated according to the prescription in (52) and compared to the threshold Fth- Cells 
in which the threshold is exceeded are marked as supercritical. 


3. 


Half the bundles in supercritical cells are moved to adjacent cells, following the direction indicated by the 
Magnus force imbalance — Fh^. 


4. The number of unpinned bundles is recorded. 

5. Steps (l)-(4) are repeated until no more supercritical cells are present. 

6. One of two things is done: in the spin down case the spin frequency of the crust (and thus |u„|) is decreased 
by bAt, with At the length of the time step, and the number of background unpinned vortices is proportion¬ 
ally reduced. Otherwise the authors mimic thermal creep by redistributing bundles in random cells to their 
neighbours. 


In general, the system requires fine-tuning of e and F^^ to obtain scale-invariant avalanches. In this sense, the 
automaton does not strictly display SOC dynamics; it does not fine-tune itself |164) . Once stationarity is achieved, 
the waiting times Tg are exponentially distributed and the distribution of sizes is indeed a power law, with exponent 
b satisfying 2.2 ^ tb ^ 5.5. Strong thermal creep leads to deviations from a power law and an excess of small events. 
In general, no correlation is observed between waiting times and the size of glitches, in line with pulsar data (see Sec. 
2). The rise times of the glitches (i.e. the duration of the avalanches) are also distributed as a power law. 

In general, the behaviour of the system reproduces the observational data well and gives an explanation for the 
collective behaviour of vortices. Furthermore the assumption that nearest neighbour interactions dominate appears 
to be confirmed, at least on a small scale, by quantum dynamical simulations of vortex motion, as we discuss in Sec. 
5.3. Naturally, the automaton above is simplified in many ways. For example the pinning threshold is taken to be 
uniform, so only an annular shell (‘critical circle’) in the star is near the unpinning threshold. Given the importance 
of this active region, a more detailed analysis should include a realistic density dependence of the pinning force, such 
as that used in the snowplow models, (see Sec. 4 for a description of how the active region is defined in these models), 
and a fuller study of how the critical circle changes radius as the pulsar ages. Further work is also needed to quantify 
the effect of the bundling parameter Bf,. As B^ increases, the size distribution develops a bump which is suggestive 
of a system-size effect. 


B. Coherent Noise and Variable-Pinning Avalanches 

The main difference between avalanches triggered by coherent-noise and nearest-neighbour processes is that the 
former occur in response to changes in the global, external driver and occur even in the absence of large-scale spatial 
inhomogeneities. The physics is as follows. If the pinning thresholds span a range, vortices unpin readily from 
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low-threshold sites and ‘stick’ in high-threshold sites over time, developing temporal correlations which lead to scale 
invariant behaviour, even when the system is spatially uniform. In the pulsar case the driver is the Magnus force that is 
built up by spin down. The pinning sites are assumed to be abundant and homogeneously distributed on macroscopic 
scales (i.e. we do not assume the existence of large scale capacitive regions with stronger pinning, as postulated in 
Sec. 5.1 and Alpar et al [251[TB5]1. Consider, for example, an ensemble of pinning sites (that are assumed to be 
defects in the nuclear lattice) each denoted by pinning strength F^, where i labels the site. The pinning strengths are 
drawn from a top-hat distribution 4>{Fp) with mean Fq and halfwidth A, such that 

= (2A)-i0(T’p -Fo + A)0(-Ap + Fq + A) (53) 


where 0 is the Heaviside step function. Each one of N pinned vortices occupies initially the site of a defect, and the 
system is evolved in time with the aid of an automaton following four simple rules m- 

1. A value of the global Magnus force Fm is drawn from a distribution '!/’(^m) which is assumed to have the same 
form as the observed glitch waiting time distribution, for a pulsar with mean glitch rate A: 

tpiFu) = exp (-Fm/ct) (54) 


with 


a = 2'KuRpnK/X 


(55) 


Fm rnay be regarded as a measure of the volume-averaged Magnus force that builds 
a more realistic model, equation (55) would involve not just A but also (1 — bvlbv^r) 
fluctuating and critical crust-superfluid lags respectively. 


up before a glitch. In 
where bv and bvcr are 


2. A small fraction / ^ 1 of vortices are allowed to unpin at random due, for example, to thermal fluctuations or 
creep. 


3. The stress Fm is applied to all the remaining (1 — f)N pinned vortices, which are allowed to unpin if F^ < Fm- 

4. Each unpinned vortex repines at a nearby defect j and is assigned a new threshold F^. 

Step (2) is necessary to avoid the system settling in an equilibrium with all vortices accumulating in ever stronger 
pinning sites [if <('(Fp) > 0 as F), —>■ oo], such that unpinning becomes rarer as time progresses. A more elaborate model 
that allows for nearest neighbour interactions would not have this requirement and could produce scale invariance 
even in the absence of thermal creep. 

Warszawski and Melatos |165j found that this model reproduces aspects of the observed glitch population. For 
example, it gives rise to a power law distribution of sizes with exponent ranging from zero to two. Fits to the seven 
most active Poisson-like glitchers suggest 0.2 < Fg/cr ^ 3 and 0.1 < A/Fq ^ 0.9. The main predictions of this model 
that are different from the vortex avalanche model are: 


• The presence of aftershocks, due to the fact that the system tends to populate sites with larger F^ as time 
passes. After a large glitch many sites are repopulated at random, leading to a lower average (Fp) in occupied 
sites, and thus to a higher probability of another glitch. 

• Glitch rise times are also distributed as a power law with exponent ~ 1. This prediction is of course very difficult 
to test, as the glitch rise time is too short to be resolved by current radio observations. 

• There is a correlation between sizes and waiting times, which is absent in the avalanche model. Such a correlation 
would be erased to some extent (to be determined by future studies) if the coherent-noise process is combined 
with nearest-neighbour interactions (Sec. 5.1) in a more realistic model as foreshadowed above. 

• Spikes may be present at the upper and lower ends of the size distribution for A ^ cr, Fg, i.e., for 4>{Fp) narrow. 
The absence of the spikes in data suggests either A « Fg or that radio timing experiments do not resolve the 
smallest glitches. Spikes also occur, in principle, if '^(Fm) has a quasiperiodic term proportional to S{Fm — Fg), 
but they do not contain much integrated probability; it is hard to extract much significance from their absence 
in the data. 
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FIG. 4: Example of results for the spin down of the crust 
in different numerical experiments, with varying ratio 
of occupied to free pinning sites, ripin/np, as described 
in Warszawski & Melatos (2011) |23| . One can clearly 
see that the spin-down is more spasmodic for a higher 
number of pinned vortices. 


C. Gross-Pitaevskii simulations 


A cellular automaton finesses the (huge) challenge of simulating the collective behaviour of the « 10^® vortices 
which thread a NS interior. However, the choice of the ‘rules’ that govern the automaton can be informed by what 
smaller-scale simulations teach us about the collective behaviour of vortices. 

To study this problem Warszawski and Melatos [IS] evolve a two-dimensional, decelerating, pinned superfluid by 
solving the two-dimensional dissipative Gross-Pitaevskii equation (GPE) 

~ ^ - (/i - (56) 

where H is the angular velocity of the frame (i.e. of the crust), is the superfluid order parameter, V{x'^) the 

external potential (which represents both the container and the pinning sites), fi the chemical potential; g parametrises 
the strength of the self interaction, and the angular momentum operator is given by Tj, = —iTizd/d(j). The term 
—^dijj/dt is a phenomenological dissipative term. 

The GPE accurately describes systems with weak interactions, such as Bose Einstein Gondensates (BECs), and 
has been successfully used for a variety of analytical and numerical studies. It is debatable whether the nuclear 
many-body forces in a neutron star are ‘weak’ in this sense; more work is needed on this question. To represent the 
typical situation in a NS (‘crust’) the potential is taken to be a trapping potential which represents the container and 
several ‘spikes’ to represent the pinning sites, of the form: 

E = Et,ap + ^ V41 + tanh(A|r - i?,|)], (57) 

i 

where Erap is the background potential, Ri is the position of the fth pinning site and E and A parametric their 
strength and width. The response of the container to vortex movement is calculated according to: 


dQ 
^ dt 


d{L,) 

dt 


+ Nem, 


(58) 


where is the moment of inertia of the crust and Nem is the external electromagnetic spin down torque. 

Warszawski and Melatos [23j carried out a series of numerical experiments with different parameters. In general, 
collective motion of vortices is observed, with vortex avalanches giving rise to glitches in the spin-down and distri¬ 
butions of sizes which obey power laws, with exponentially distributed waiting times. The main conclusions can be 
summarised as follows: 


1. Vortices unpin in ‘herds’ under the action of two knock-on processes: the proximity effect (unpinning due to 
vortex-vortex interactions) and acoustic knock-on (unpinning to sound waves shaking the vortex). The latter 
dominates the former when 7 is small and biases the dynamics towards the coherent noise model (Sec. 5.2), 
whereas the proximity effect dominates when 7 is large and the dynamics is more SOG-like (Sec. 5.3). This 
validates earlier suggestions by Alpar [2511168] that random unpinning in capacitive vortex regions of the crust 
could initiate a vortex avalanche. 










Models of Pulsar Glitches 


20 




FIG. 5: Probability density function of fractional glitch sizes (left) and cumulative probability function for waiting times 
between glitches (right), obtained from the simulations of Warszawski & Melatos (2011) |23| with varying ratio of occupied to 
free pinning sites, ripin/np shown in figure]^ The grey line on the left hand graph is a power law least square fit to the solid 
histogram, with exponent -0.801. 


2. Vortices tend to skip over pinning sites and more radially outwards by a few Feynman distances before knocking 
on a neighbour, thereby maintaining a roughly (but not exactly) regular Abrikosov lattice. The pinning strength 
plays a fundamental role, with stronger pinning giving rise to larger glitches, and a broader range of potentials 
giving rise to a broader range of sizes and more frequent glitches |2S] 

3. A strong external spin-down torque leads to smaller, more frequent, glitches than a weak torque 

4. A heavier crust produces more large glitches and longer waiting times than a lighter crust 

5. Waiting times are longer and glitches larger when vortices and pinning sites are equally abundant, than if pinning 
sites are more abundant than vortices. Further increases in pinning site abundance do not change the results, 
as vortices tend to simply move across many pinning sites to readjust the Magnus force 

In figure]^ we show an example of different spin down curves for the crust obtained with different setups of the GPE 
code, with varying ratio of pinned to free vortices, npin/np, as described by Warszawski & Melatos [23]. In figure 
l^we also show the probability density function of fractional glitch sizes and the cumulative probability function for 
waiting times between glitches, obtained from these simulations. 

Motivated by the GPE results, Warszawski & Melatos |21| developed an asynchronous automaton (i.e. an automaton 
in which the time step is itself a stochastic variable that depends on the state of the system, see eg. |169j for 
a description of how this choice influences the system) to study the statistics of pulsar glitches. They show that 
a simple Poisson process in which vortices unpin independently (e.g. by thermal creep) cannot explain the broad 
distribution of glitch sizes observed in the pulsar population, even if one accounts for variable pinning strengths (with 
the only difference in this case being that stronger pinning sites are preferentially occupied, leading to occasional 
large glitches). In order to obtain a broad power-law size distribution knock-on effects need to be included (i.e. one 
has to allow for the possibility that upon unpinning a vortex will initiate an avalanche). Nevertheless also in this 
case a power law distribution is only obtained if the spin-down rate and unpinning rate (which can be expected to 
be a function of temperature if we allow mainly for thermal creep) are fine tuned. This fine tuning is unlikely in an 
astronomical setting and further developments are required to understand this issue. 

Finally let us note that the main limitation in applying the results discussed above directly to an astrophysical 
setting is linked with the scale of the simulations. In general in all the examples above vortices are separated by a 
few pinning sites, and the overall superfluid flow is strongly influenced by motion of the vortices themselves (hence 
the importance of knock-on effects). In a realistic neutron star the inter-vortex separation is much larger than the 
average distance between pinning sites, with an average vortex encountering up to « 10^° pinning sites between its 
equilibrium position and its nearest neighbour. Two things can then happen. On the one hand, an unpinned vortex 
may skip over « pinning sites before repinning. This is not as unlikely as it sounds, as the systems constantly self¬ 
adjusts to keep the effective Magnus-modified pinning potential marginally stable and roughly uniform everywhere. 
Alternatively vortices could be too sparse for proximity effects to play an important role and lead to knock-on effects. 
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In this case vortices will creep out gradually, as we describe in the following section. The exact nature of the motion 
of a realistic vortex in a realistic neutron star thus remains an open, and central, problem. 


VI. HYDRODYNAMICAL INSTABILITES 


An important ingredient of most of the models described above is a pinned superfluid in the crust or core which 
can store angular momentum by lagging behind the ‘normal’ fluid and the crust. An important implication of this 
picture, which has received little attention until recently, is that the superfluid analogue of the two-stream instability 
can occur imisH], possibly triggering vortex unpinning and hence a glitch. 

One mode of oscillation triggered this way is the r-mode, a toroidal oscillation where the restoring force is the 
Coriolis force. This particular class of inertial modes has already attracted much interest as it can be driven unstable 
by gravitational wave emission |170L 1171) . Its velocity field is particularly simple in the two-fluid description 11721 , 
even when one accounts for a lag [2^1173) . In particular we study linear perturbations of the equations in (14)-(16l, 
with a rotational lag A = (Up — fln)/np in the background, as would be the case if vortices are pinned (B = 0 and 
B = 1 for strong pinning). One can solve for a velocity field of the form (where x = p, n labels the fluid) |173) : 


SvL = -- 


im 


’ sind 


UlYrel + 


* sin 9 


ULdeVre^. 


(59) 


where U^{r) are the mode amplitudes and Y^{9^(j)) are the standard spherical harmonics. Consider the simple 
case of vanishing entrainment and a constant density star, for which the equation of state of the perturbations is 
incompressible, i.e. (Vidv^ = 0). One can find simple solutions for the I = m case, for which and the 

mode frequency a takes the form [173] : 
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a = — 


(m -f l)a;, 


-[1-Xp + A±D^^% 


D = (1-I-Xp)^-I-2A{1-I-a;p[3 - m(m-I-1)]} 


r-modes are therefore unstable for 


The growth time for the instability is 


m > rric = 
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\/2xpA 
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(60) 

(61) 

(62) 

(63) 


where P is the rotation period of the star. The growth time scale must be compared to the time scale on which shear 
viscosity damps the mode, which, for a constant density star with i? = 10 km and M = 1.4M0, is given by m 
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that is, the mode grows provided that 
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This constraint, together with (62) tells us that the mode first becomes unstable at the critical lag 
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(64) 


(65) 


( 66 ) 


Glampedakis and Andersson |29j used angular momentum conservation to compare (66) to maximum glitch sizes and 


find Ac < (/p//n)(A/flp), suggesting that this mechanism can set the maximum size of a glitch, and larger lags are 
prevented by the instability (which may, in fact, trigger the glitch). 

In figure ([^ we show a comparison of the prediction for the critical lag with data, for an assumed value of the ratio 
Ip/In ~ 0.02. The model can explain the different maximum sizes in different pulsars by accounting for differences 
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FIG. 6: A comparison of the critical lag 
in (661 against data for glitches of several 
pulsars, with a given spin period P. The 
red dots are systems where the tempera¬ 
ture can be measured, while for the trian¬ 
gles the temperature is estimated from a 
simple modified URCA model to be T « 
3.3 X 10®(ts/lyr)“^/® K [^. The unsta¬ 
ble region is in gray and rric represents the 
critical m for the onset of the instability. 
For higher m there is a range of unsta¬ 
ble inertial modes. A ratio Ip/P ~ 0.02 
is assumed and the ratio Aflc/flc used in 
Glampedakis and Andersson is equivalent 
to A/f2. 


in temperature (which is inferred from standard cooling models |29| 1 and age. A more complex analysis, including 
entrainment m, reveals that not only does the instability persist, but that a secular (long timescale) instability 
exists for low I in the pinned regime. 

Link has also shown, in a plane wave analysis, that instabilities still exist, both in the core and the crust, if one 
does not assume perfect pinning but allows vortices to creep out, i.e. a regime such that B << S [m 1175) . The 
instability is thus unlikely to be stabilised by viscosity or mutual friction and can play a significant role in glitches. 

Another important instability in laboratory superfluids is the Donnelly-Glaberson instability |176) . In superfluid 
helium this signals the transition to turbulence. When counterflow along the rotation axis exceeds the threshold 
wdg = 2-\/2^r2, Kelvin waves destabilise the vortex array and create a vortex tangle. For standard parameters 
\wdg\ ~ 1.24(0.1 s/p) ' , where v = («;/47r) log6o/ao, with 6o the inter vortex spacing and oq the radius of the 
vortex core pH) . 

In a realistic neutron star the crust is spun down first by electromagnetic torques, and the rest of the fluid follows 
in response to magnetic torques or Ekman pumping. Peralta et al. [281 1177) showed that, in the absence of magnetic 
fields, a sizeable axial counterflow develops which exceeds the threshold for the Donnelly-Glaberson instability in most 
of the star. The mutual friction force changes form and strength in those regions where the threshold is exceeded. 
Note that these calculations are carried out in the Hall-Vinen-Bekarevich-Khalatnikov (HVBK) [17811179) formalism 
that describes superfluid helium, i.e. a superfluid condensate and its massless excitations. Nevertheless the results 
should not depend strongly on this choice and highlight the possibility that turbulence is likely to play an important 
role in pulsar glitches, with transitions between laminar and turbulent mutual friction consistent with the short spin- 
up and long inter-glitch timescales observed [55]. Recently Melatos and Link [180) have also shown that shear-driven 
turbulence exerts a fluctuating torque on the crust, which may contribute to the timing noise observed in pulsars. 

Mastrano and Melatos [181) suggested that a Kelvin-Helmholtz instability at the interface between the and ^P 2 
superfluids in the star may trigger the glitch, complementing the bulk two-stream instability discussed above. 


Relaxation 

The long timescales (weeks to months as described in section 0 observed in post-glitch relaxation were the first 
evidence that neutron stars contain a superfluid component, loosely coupled to the crust [18]. Two classes of models 
have been advanced. Alpar and collaborators [2111182) suggested that, by analogy with flux creep in superconductors, 
vortices ‘creep’ out as the star spins down, and the response of the creep rate to a glitch determines the spin rate 
of the star and gives rise to the observable relaxation. Alternatively, if pinning is weak in the crust, as suggested 
by Jones [78) . the vortices move with the neutron superfluid, and mutual friction is the dominant mechanism in the 
relaxation . Van Eysden and Melatos [85| analysed the general case where the ratio of the Ekman and mutual friction 
time-scales can be small or large. They obtained good agreement with Grab and Vela data as well as laboratory 
experiments with superlfuid helium (to an accuracy of 0.5 % in the latter case [116) 1. From the analytic solution of 
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FIG. 7: The inferred superfluidity coefficients BE 
(where E the Ekman number, i.e. E = 
with rj the shear viscosity coefficient) plotted against 
Pn(l + K), where K is the ratio of the total moment 
of inertia of the fluid to that of the crust (See van 
Eysden & Melatos [ 85 ] for details). The lightly and 
dark shaded regions correspond to the range inferred 
from Vela and Crab data, respectively. The theoretical 
curves are calculated by van Eysden & Melatos [ 85 ] for 
three different temperatures: T = 10^ ® K (top, red), 
T = 10 ^ K (centre, green) and T = 10 ®'® K (bottom, 
blue). In these models Pn is varied for a fixed value of 
K = 50 . 


van Eysden and Melatos [85] one can also obtain a general condition for observing an ’overshoot’ in the relaxation, 
i.e. for observing the post-glitch frequency to drop below the value it would have had in the absence of a glitch and 
subsequent relaxation, as is the case for the Crab, but not Vela. In the model of van Eysden and Melatos [SS] this is 
found to occur if immediately after the glitch one has flp < fin < + AOgi be. after a glitch of size AEIg the crust 

is spinning slower than the superfluid component. The different behaviour of the Crab and Vela pulsars would then 
suggest a difference in the internal vortex and/or pinning site distribution. The viscosity coefficients of nuclear matter 
obtained by fitting the relaxation of both pulsars are, however, consistent with each other, and with the theoretical 
predictions of van Eysden and Melatos [85] , as can be seen from figure 


VII. VORTEX CREEP 

In a seminal work Alpar et al. |2Ij . inspired by the analogy with flux creep in type II superconductors, developed 
a theory of non-linear coupling between the crust and the neutron fluid, due to vortex ‘creep’ between pinning sites 
and applied it to the recovery of Vela giant glitches |182j . The main idea is that the core superfluid couples rapidly 
to the normal component, due to electron scattering on vortex cores [S^, according to The observed dynamics 
must then arise from the crust. This conclusion was supported by the observed increase in spin-down rate after the 
glitch, which points to a few percent of the moment of inertia being decoupled, the appropriate amount for the crust 
superfluid. Short term increases in the spin-down rate observed in later glitches, however, can exceed 100%, implying 
that ‘decoupling’ the superfluid in the crust is not enough; one needs a time-dependent internal torque (e.g. a change 
in mutual friction mechanisms) |11511183] . 

Vortices in the crust pin to ions [IlEl]- On long time-scales, however, the neutrons follow the spin-down of the 
crust, as vortices migrate from one pinning site to the next, via thermal excitations or quantum tunnelling |184j . i.e. 
they ‘creep’ out. The Magnus force biases the motion in the outward direction, leading to an average radial velocity 
(vr) 7 ^ 0. Suppose one views the Magnus force as arising from a fictitious pinning energy [?!] 

where Aflc is the critical lag for unpinning (i.e. after which the pinning force can no longer balance the Magnus 
force). Eor outward radial motion, the energy barrier that a vortex must overcome is lowered to 

E+= Ep-AEp = Ep{l-An/A^a), (68) 

while for motion in the —r the barrier is higher: 

E_ = Ep + AEp = Ep{l + An/An^). (69) 

In other words for a vortex to move radially in work has to be done against the sum of the pinning force and the 
Magnus force that tends to push the vortex out. For a vortex to move out, on the other hand, work has to be done 
against a lower total force, given by the pinning force minus the Magnus force. 
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In thermal equilibrium the radial bias in the pinning potential leads to more vortices randomly unpinning and 
moving out than in, and thus to an average outward creep speed: 


Vr 


Vo 


exp 




Vo exp 



(70) 


where vg is the typical velocity of microscopic motion of the vortex lines between pinning centres, which Alpar et 
al. [21] take to be due to the Bernoulli force, with vg « 10^ cm/s. Note that this assumed value is critical to the 
calculation, and is simply taken to be a constant in the original model. Several authors have estimated the average 
speed of free vortices, and found it to vary significantly in different regions, and depart from the constant value given 
above (HU [M] 1185] , leading to significant differences in the predictions of the theory, as we shall see in the following. 
More detailed calculations of E± |186j account also for the tension of the vortex line; vortices behave rigidly in weaker 
pinning regions, unpinning simultaneously from several sites, while they are more flexible in stronger pinning regions 
(in which the pinning force dominates over the tension), with unpinning dominated by single-site events. The relative 
importance of the tension can be expressed in terms of a ‘stiffness’ parameter Ts = Trg/Fpl, where T denotes the 
tension, Fp is the maximum pinning force, rg is the range of the pinning potential and I is the internuclear spacing [186]: 

Pn \ f Ep / rp \2 / I 

lO^^^g/cm^y ylMeV/ VlOfm/ \50fm 

In the limiting cases of AH « Aflc and Ail << Ailc, the activation energy Fa (not to be confused with the pinning 
energy per site Fp) to unpin can be written simply as: 




Fa - 


for 

Ail « Ailc, 

(72) 

Fa - 

» b.lFp^/T~s - 

Ail ' 
Ailo^ 

j for AD « Ailc 

(73) 


In order to obtain a radial creep rate Link et al. |184j also estimate the rate at which vortices ‘attack’ (either due 
to thermal excitations or to quantum mechanical zero point motion) the barrier and subsequently repin. They find 
that 


Vr 

Vr 


UJs . /'2TTkTeff 


27r 


o-KAFss 


1/2 


for Tg < 1 


for Tg > 1 


(74) 

(75) 


where the asterisks indicate that the quantities are evaluated for the characteristic number of pinning bonds j* 
involved in an unpinning event, Wg and ujf are the characteristic Kelvin oscillation frequencies for stiff and flexible 
vortices I min is the minimum distance between pinning sites, and one has T^ff = TqCoih{Tq/T), where Tq is 

the crossover temperature above which vortices surmount the barrier thermally, and below which quantum tunnelling 
through the barrier dominates. 

In order to make quantitative comparisons with glitch recoveries Alpar et al. |182j assume that there exist physically 
distinct regions of different pinning strengths. Subsequent calculations do, in fact, show variations of pinning energies 
with density, and random variations of the pinning force due to different orientations of vortices [SDJ I15HL I187j . 
Depending on the temperature, a given region with a given pinning strength, can be in either the linear or non-linear 
regime of vortex creep |188j . If Fp exceeds the binding energy of a nucleus to its equilibrium site, so that it is 
energetically favourable for a vortex line to dislodge nuclei, one has ‘strong’ pinning. When this is not the case, but 
the vortex core radius ^ (the superfluid coherence length) is smaller than the lattice spacing, one has ‘weak’ pinning. 
Finally, when the superfluid gap decreases and the vortex core is large enough to encompass several nuclei, there is 
little change in energy between nearby configurations (as the number of pinned nuclei remains approximately constant 
as the vortex moves) and one has ‘super weak’ pinning. 

The equations of motion for the system are 


7pDp — 


No: 




(76) 


where ilp{t) is the angular velocity of the crust and all components tightly coupled to it, iln{t) is the angular vel ocity 
of the superfluid in the pinning region, and Vr is obtained in the pinning regions from the expression in (701, 
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with An* = flp — The external spin-down torque is given by Next and /„ is the moment of inertia of the crust 
and all components tightly coupled to it, which in the formulation of Alpar et al. [2T] includes the entirety of the core 
neutrons, which are assumed to be coupled on timescales much shorter than those associated with the glitch [53]. The 
glitch trigger is not specified; it is introduced as an initial condition. Elsewhere Alpar et al. [3^ discuss how random 
unpinning of vortices in ‘capacitive regions’ (i.e. stronger pinning regions with an overdensity of vortices) can lead to 
vortex avalanches, a suggestion that has indeed been verified by GPE simulations [24] . 

Given an initial angular velocity perturbation Ar2i(0) induced by a glitch there are two regimes of dynamical 
response of the equation in (76). If the temperature is sufficiently large compared to the pinning energy, then a steady 


state (in which the neutrons spin down at the same rate as the crust) can be achieved with a small enough lag that 
the response of the spin-down rate to the glitch is linear in the size of the perturbation. The contribution of the ith 
pinning region shows up in the post-glitch response AQp as an exponential decay: 


= - 


A Af2,(0)e 




ns 


(77) 


where li is the moment of inertia of the region, / is the total moment of inertia of the star, and the local coupling 
timescale in the is given by: 


n,i = 


kT 

a; 


(78) 


If, on the other hand, the temperature is low compared to the pinning energy, the steady state lag is large and the 
response to the perturbation is non linear. In this regime the effect of a region k on the response of the spin down 
rate has a Fermi function dependence on the glitch size, of the form: 


Afl 


p,fc 


it) = 


1 - 


1 -I- [exp{to,k/Tn,k) - 1] exp(-t/T„,fc) 


(79) 


where Ik is the moment of inertia of the region, t^^k 


Ar2(0)/|I7n| and the relaxation time can be expressed as: 


kT Ane 
Ep |J7n| 


(80) 


Greep islinear if ti < with the transition occurring at a pinning energy of Ep/kT « 31 [188] . corresponding to 
timescales r « 1000 days for Vela, for weak pinning. Given that this is roughly the periodicity of Vela glitches it 
is clear that non-linear relaxation will be difficult to observe. However, as we shall discuss in the following section, 
the timescale is much shorter if the transition occurs in a super weak pinning zone, and potentially observable. In 
fact most of the timescales involved in the relaxation of Vela glitches are interpreted as linear response of different 
regions, with only the longest timescales open to the possibility of a non-linear response (although combinations of 
linear responses are possible) |189| . The observation of a Fermi-function like behaviour would, on the other hand, 
suggest the presence of non-linear creep in the star. 

To explain the observations of the Vela pulsar Alpar and collaborators develop a ’minimalist’ glitch model [21111881 
m\- They assume that vortices accumulate in a ’capacitive’ region Al, which is adjacent to a depletion region B, in 
which no vortex motion is present (i.e. is not contributing to the long-term spin down). During the glitch vortices 
unpin from region Al, pass through region B and finally repin in a region further out, A2. These are the ’active’ 
regions, while there are ’passive’ regions further in the star in which no vortex motion occurs and that react to the 
glitch only via the change in lag induced by the glitch (i.e. for these regions AD = ADgiitch)- 

After integrating the contributions for these different regions, one obtains the general form for the post glitch 
response as [189] : 


^ = (t + t) - (t + t) + 

^ r dh AD;(0)e-*/"‘ ^ /■ ^ A_ I _ 

J I IDpT/ J I \ 1-b [exp(to(t)/T„ -1] exp(-t/T„) 

where the subscript I refers to regions in the linear coupling regime and nl to regions in the nonlinear coupling regime. 
IA is the moment of inertia of region Al (the region from which vortices depin, Ib of the depletion region through 
which vortices travel, while Ia 2 « Iai and is neglected. One has tg = AD(0)/|Dn| for regions A and B. 

By fitting the expression in Q to the post-glitch relaxation of eight Vela glitches, Alpar et al. [189] reach the 
following general conclusions: 
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• Short timescale exponential relaxation generally corresponds to the response of passive regions (i.e. regions 
through which no vortex motion occurs) in the linear regime, which are identified with super weak pinning 
regions with pinning energies of around Ep « 0.3 MeV. 

• The full non-linear Fermi function response is not observed in the recoveries of Vela glitches 


• The remaining longer timescales and linear terms in the relaxation are more difficult to interpret, given that 
for to << T the nonlinear term reduce to an exponential. If one has relaxation timescales r of approximately a 
month and glitches of approximately Aflp « 10“^ rad/s, then to ~ lOdays << r and thus we are in the regime in 
which the nonlinear response will give rise to an exponential. One could then be observing the linear recoupling 
of two different regions, one with r « 30 days, and the other with r > 1000 days. Such a long timescale would 
lead to what appears to be a permanent offset on shorter timescales. However Alpar and collaborators |189j 
suggest that in reality one is observing non-linear recouping in super weak pinning, with regions through which 
no vortex motion has occurred giving rise to the 30 day exponential relaxation and regions through which vortex 
motion occurred giving rise to the permanent increase in spin down rate. This places the unpinning regions 
close to the transition between linear and non-linear response, and the transition between super weak and weak 
pinning at slightly lower densities. 

• Alpar et al. |168L I190j also suggest that the different glitching activity of the Crab is due to the star not 
having yet formed a fully developed network of capacitor regions, and that crust quakes may be the trigger for 
unpinning for this pulsar (as a opposed to vortex avalanches for the Vela pulsar). 


In general the theory is successful in describing the relaxation of Vela glitches with consistent parameters. Although 
there are some difficulties in explaining large increases in the spin down rate on short timescales of a minute or so 
after the 2000 and 2004 Vela glitches m, we will see in the following that these can be overcome if one allows 
for part of the core to decouple after a glitch. This is also necessary, as described above, if one is to explain the 
activity of the Vela pulsar while accounting for strong entrainment in the crust. Another important point is that 
most of the recovery can be described by exponentials, i.e. by creep in the linear regime. In this regime creep can be 
approximated in the two-fluid picture in terms of a fraction of free vortices |185j . allowing us to apply the formalism 
described in section (IB). Nevertheless the presence of several free parameters to account for different pinning strength 
and coupling timescales leads to a theory that is, in fact, difficult to falsify. Recently Link m has examined the 
problem of vortex creep in more detail His conclusion is that both the linear regime and super weak pinning invoked 
by standard vortex creep theory described above are unlikely to be realised in practice. Let us review Link’s revised 
theory of vortex creep, dubbed vortex ’slippage’. 


A. Vortex slippage 


Recently Link [ST] has revised the vortex creep model, by generalising two of the main simplifying assumptions of 
the original theory. The first is the form of the pinning barrier, for which Link |81j uses more realistic values, such as 
those in equations (73), the second, is the typical velocity of an unpinned vortex, which in the original theory is fixed 
at Vq ~ 10^ cm/s. In vortex slippage theory the velocity of a free vortex moving from one pinning site to the next is 
obtained by accounting for tension and drag forces, resulting in an average velocity for a vortex of the form [81] : 


< V >= rAH ( ^ sin 20(^7 -|- cos^ 9d4> I e 


-{EJkT} 


(82) 


where dd is the dissipation angle, related to the mutual friction coefficient by 2B = sin 20^6 The post glitch 

response is found to depend non-linearly on the size of the glitch, with pinned regions responding on a time-scale 
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where Gd is the spin down timescale and AHg the size of the glitch. Exponential decay can only occur on timescales 
longer than td, and thus the shorter exponential timescales observed in the recovery of Vela glitches are unlikely to 
be due to vortex slippage in either the crust or core. Haskell et al. I22| interpret the latter timescales as due to free 


vortices coupled by mutual friction in the outer core (see section VHI). Finally Link estimates the pinning strength, 
also accounting phenomenologically for strong entrainment, and the possibility of pinning both to the crustal lattice 
and to flux tubes in the core. His calculations, together with the estimates for the coupling timescale in (83) suggest 


that the linear creep regime introduced by Alpar and collaborators pTl 118811189] is not realised in practice and that 
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super-weak pinning is is precluded by thermal fluctuations. Observations of a delayed response in the recovery after a 
glitch, over a timescale proportional to the glitch size as in (83) would, on the other hand, would provide a hint that 
vortex slippage is occurring in neutron stars. 

Finally let us note that vortex creep is a dissipative process and creep heating ll91H193j might explain the 
unusually high (for their age) temperatures that have been measured for some old, recycled millisecond pulsars, such 
as J0437-4715 [l94] . 


B. Thermal glitches 


Temperature clearly plays an important role in the creep model. In the previous discussion however the trigger 
mechanism for the glitch itself is left unspecified and one does not account for heat deposition during the glitch and 
how this would affect the motion of vortices. This problem was examined by Link and Epstein |195j who considered 
the effect of depositing heat in a region where the glitch is triggered (e.g. due to the strain released by a quake, 
or dissipation due to mutual friction during a vortex avalanche). If the trigger is a crust quake, one could deposit 
up to 10^^(dc/10“^) erg in the crust |195) . After the trigger deposits heat locally, a heat front travels through the 
crust, increasing the creep rate and spining-up of the crust further. The glitch then ends either when the thermal 
energy diffuses or when the superfluid and crust velocities are equal. This scenario was further considered by Larson 
and Link |124j who simulated Crab and Vela glitches by solving the equations of motion for creep in (76), with the 
detailed prescription for the creep velocity in (74]|75), coupled to the thermal diffusion equation 


dT 

Cv-^ = Vj(KTV*r) + hf - 


(84) 


where c„ is the specific heat, kt the thermal conductivity, tu the neutrino emissivity and hf the heating rate per unit 
volume due to superfluid friction in the crust, which corresponds simply to the rotational energy dissipated in the 
crust due to vortex motion mmss] and is obtained from H — J c?/nAn|lin|, dividing by the volume of the crust. 
The simulations are spherically symmetric and initiated by depositing heat (corresponding to typical energies of 10^^ 
erg) in a spherical shell centred on a density of p = 1.5 x 10^^ g/cm^ with a width of 40 m. 

The results of the simulations are compared to glitches in several pulsars, e.g. Crab and Vela. The agreement is 
good, with the interesting conclusion that, for similar heat depositions, glitches in hotter pulsars (such as Crab) are 
generally smaller and slower than in older pulsars (such as Vela), essentially due to the larger specific heat |124j . This 
model thus accounts for the continued rise observed in the 1989 Crab glitch for a day (but is not commonly observed), 
which the simulations of Larson and Link do not reproduce if the glitch is introduced as a sudden jump in frequency 
with no heat deposition. Furthermore the large heat deposition can lead to thermal pulses at the surface, which are 
close to the current observational upper limits on post-glitch heating of Vela (see section E- 


VIII. MUTUAL FRICTION 


The creep paradigm relies on vortices being pinned in the crust, and only a small number ’creeping’ out. Clearly if 
the pinning is very weak and most vortices are free this will not be the case. This was suggested by Jones [TS], who 
pointed out that for a infinitely long vortex the pinning interactions from the surrounding ions would average out, 
leading to it being essentially free. He asuggested that this also happens for a realistic, rigid vortex, and that vortices 
in the crust are close to co-rotation with the neutron superfluid. The observed recovery timescales are then due to 
mutual friction. Although recent calculations reveal that, for realistic tension forces, there is still a sizeable pinning 
force [SD], the fact that creep is expected to be in the linear regime for most of the post-glitch relaxation lets one 
simulate a glitch using of the multi-fluid formalism described in section (IB). Haskell et al. [12] have integrated the 


equations in (14 16) in cylindrical symmetry by averaging the mutual friction coefficients over the length of a vortex. 
Mutual friction is t ake n to be due to electron scattering in the core, with a mutual friction coefficient of the order 
B « 10“^ given by (11), and to phonons in the crust, with a mutual friction coefficient of the order B « 10“® [Ml 165] . 

The system is evolved with both fluids originally co-rotating and all vortices in the crust pinned, i.e. B = 0 (note 
that Haskell et al. |22| ignore the effect of entrainment, which reduces the lag accumulated during spin down). As the 


evolution proceeds the lag AH is compared to a realistic pinning profile (such as those described in section IV) and 
if it exceeds the critical lag for unpinning AHc mutual friction is switched back on, mimicking the effect of vortices 
unpinning. The region then relaxes back to equilibrium and no re-pinning is allowed until after a glitch. Once the lag 
reaches the maximum of the critical lag (we remind the reader that the critical lag varies with density, and thus has a 
maximum in the crust, as described in section IV), a glitch is initiated. This is done by assuming that as vortices move 


out they excite Kelvin waves, thus leading to strong dissipation and rapid recoupling. In practice this is achieved by 
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FIG. 8: LEFT: Comparison of the short-timescale post glitch relaxation simulated by Haskell et al. |22| for varying values of 
the drag parameter TZ in the core, to the observed relaxation of the Vela 2000 and 2004 glitches. The pre glitch spin down 
has been subtracted. In general we can see that larger values of TZ are favoured, with TZ ~ 10“'*, as expected from electron 
scattering of vortex cores, consistent with the data. Values TZ < 10“®, as would be the case if a large number of vortices remain 
pinned to flux tubes in the core, are not consistent with the data, given this model. RIGHT: Example, in the setup of Haskell 
and Antonopoulou |119| of the different relaxation behaviours one can have by increasing B during the glitch, in an 80 m region 
of the crust, with a background (pre-glitch) S = 5 x 10“®. As B increases (and the rise time decreases) passes from 0.9% 

for B = 8 X 10“^, to 1.3% for B — 9.5 x 10“^, and finally for S = 2 x 10~® an exponential trend begins to appear. 



FIG. 9: The left panel shows a schematic representation of the coupling in the crust and outer core. In the middle panel we 
show the large decrease in coupling timescale that would lead to a glitch that decouples part of the core, and then relaxes (such 
as a giant glitch in the Vela). In the right panel we show the situation for an avalanche that frees only a small number of 
vortices and thus only involves the crust. In the picture of Haskell and Antonopoulou m this gives rise to glitches for which 
the relaxation appears essentially as a step in frequency and frequency derivative, with no exponential relaxation. 


switching the mutual friction parameter to Bk ~ 10“^ in the region near the maximum that hat has not yet relaxed to 
its equilibrium lag, i.e. for which AQ > VL-p/2VLaB, where B is the crustal mutual friction parameter before the switch 
to Bk- This leads to a rapid rise (typically of the order of seconds) and glitches of the order of Afl/n « 10“® — 10“^. 

In hgure we show fits to the Vela 2000 and 2004 glitch recoveries. In general the results are consistent with 
electron scattering in the core as the main source of mutual friction, but not with the drag being much weaker, as 
could be the case if vortices pin to flux tubes in the core during the glitch. 
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Haskell and Antonopoulou |119) adapted the above model to smaller glitches, to show that random unpinning will 
lead not only to glitches of different sizes, but also to glitches with different rise times. This follows from introducing 
a fraction of unpinned vortices 7 in the equations of motion ( 14pT ) as an effective mutual friction coefficient B = jB. 
The coupling timescale between the fluids is thus a function of 7 and is given by t « 1/(2HH). Large values of 7 
(i.e the unpinning of many vortices) produce preferentially larger glitches that decouple part of the core, as shown 
in figure (and as in the previously described simulations of large Vela glitches), and are followed by exponential 
relaxations. Small values of 7 lead to glitches that do not involve the core and look like simple steps in frequency 
and frequency derivative, as can be seen in figure This picture is also applied to relaxations of magnetar glitches, 
for which large increases in the spin-down rate can be obtained on long timescales simply because the magnetars are 
much slower than ordinary pulsars, leading to the recouping timescales for the interior regions r « 1/(2VLB), which 
give rise to the recovery, being consequently increased. 

An effect that is ignored in the above picture is that of Ekman pumping, which will act to communicate changes 
in the spin of the crust to the bulk fluid of the interior, van Eysden and Melatos showed that the timescales can 
be comparable to the shortest in glitch relaxations in the case of non-magnetized stars [85], but much shorter in the 
presence of magnetic fields jSS] . In the non-magnetised case the problem can be solved exactly in different geometries 
for varying ratios of the Mutual Eriction and Ekman timescale |197j . However, in the magnetic case, the dynamical 
response of the fluid is complex, exhibiting various modes of oscillation after a spin up (down) event, which may be 
observable in post-glitch relaxations [ 86 ] . 


IX. GRAVITATIONAL WAVES 

Gravitational radiation is neither scattered nor absorbed. It is therefore an ideal probe of dynamical processes 
in neutron star interiors, which often lack direct electromagnetic signatures. The next generation of ground-based, 
long-baseline, gravitational-wave interferometers are on schedule to commence science operations in the next few 
years m, with sensitivites ~10 times greater than previously achieved. It is therefore timely to review briefly the 
possibilities for multi-messenger (e.g. simultaneous radio and gravitational-wave) studies of pulsar glitches in the 
coming Advanced Detector Era. 

Gravitational wave signals from a glitch divide into two types: broadband bursts from the spin-up event itself, 
and narrowband continuous-wave emission from the inter-glitch interval and post-glitch recovery. The burst signal is 
triggered by whatever non-axisymmetric dislocation triggers the spin-up event, e.g. crust cracking or superfluid vortex 
unpinning. The energy released can excite the acoustic and intertial stellar oscillations, e.g. f and p modes [19911201] 
and their superfluid versions generating an oscillating mass quadrupole moment lasting < 0.1 s. The energy 

conversion efficiency is poorly known [ 20511205 ] . Alternatively, the rapid motion of the unpinned vortices, before 
they repin, can produce a time-varying current quadrupole moment, again on the spin-up time-scale [205]. Gurrent 
quadrupole radiation is also emitted by glitch-triggered r-modes [ 200 ] . The continuous-wave signal is generated by 
residual non-axisymmetries in the superfluid velocity field, some parts of which are erased on the post-glitch recovery 
time-scale. Elow non-axisymmetries arise from Ekman circulation [20711208] and/or oscillatory vortex structures at 
high Reynolds number, e.g. Taylor-Gortler modes [2811209] . The emission is through the current quadrupole channel. 
However, if glitches arise from vortex avalanches, then some part of the above non-axisymmetries persist beyond the 
recovery stage and the pinned vortex distribution following a glitch is non uniform, as in other SOC system [138] . A 
general argument based on angular momentum conservation in the context of glitches triggered by vortex avalanches 
shows that the persistent, current-quadrupole wave strain Hq from a glitching pulsar satisfies [ 210 ] 

/lo > 3 X 10“^^(D/10^ rad s“^)^((i/lkpc)“^(/p//„)(AD„aa:/f^), (85) 

where d is the distance to the source, and Ailmax/^ is the maximum fractional glitch size over the current spin-down 
time-scale, whether or not that largest glitch occurred during the gravitational wave observation. 

At the time of writing, the most sensitive searches for gravitational radiation from glitches have been performed 
with Science Run 5 data from the Laser Interferometer Gravitational wave Observatory (LIGO). A Bayesian analysis 
of the August 2006 glitch in the Vela pulsar, incorporating prior information about known instrumental artefacts, 
placed a 90% confidence upper limit of /iq < 1.4 x 10“^® on the wave strain amplitude of quasinormal mode ring-down 
signals with frequency 1-3 kHz and decay time 0.05-0.5 s [21111212] . The same search, conducted with Advanced 
LIGO, would be ~ 10 times more sensitive. Advanced detector searches for post-glitch, continuous-wave signals are 
also planned using algorithms for long-duration transients with phase drifts, e.g. Bayesian window template [213] 
and frequency-time cross-power EH]. If the quadrupole is generated by pinned vortices, it co-rotates with the crust, 
hence is phase locked to the radio ephemeris [ 210 ] . 

Finally, looking even further into the future beyond the Advanced Detection Era, instruments like the Einstein 
Telescope are being designed with enhanced Newtonian noise suppression in the 10-50 Hz band, where many glitching 
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pulsars reside. The prospects for further enlargement of the glitch discovery space are therefore bright. 
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